diff --git a/.Rbuildignore b/.Rbuildignore
index 7934390..3ed6d02 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -14,3 +14,4 @@
 ^src/\.clang-format$
 ^src/core/\.clang-format$
 ^src/core/test$
+^.vscode$
\ No newline at end of file
diff --git a/NAMESPACE b/NAMESPACE
index 84f46d6..f29e724 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method("*",DiffOpObject)
+S3method("*",DifferentialOp)
 S3method("*",Function)
 S3method("*",FunctionGrad)
-S3method("+",DiffOpObject)
-S3method("-",DiffOpObject)
+S3method("+",DifferentialOp)
+S3method("-",DifferentialOp)
 S3method("st_crs<-",Domain)
 S3method(plot,Function)
 S3method(plot,Mesh)
@@ -25,21 +25,10 @@ export(dot)
 export(grad)
 export(laplace)
 export(read_mesh)
-exportClasses(BasisFunction)
-exportClasses(DiffOpObject)
-exportClasses(DiffusionOperator)
-exportClasses(Domain)
-exportClasses(Function)
-exportClasses(FunctionGrad)
-exportClasses(FunctionSpace)
-exportClasses(Mesh)
-exportClasses(Pde)
-exportClasses(ReactionOperator)
-exportClasses(TimeDerivative)
-exportClasses(TransportOperator)
 exportMethods(contour)
 exportMethods(dt)
 import(Matrix)
+import(R6)
 import(Rcpp)
 import(RcppEigen)
 import(methods)
diff --git a/R/domain.R b/R/domain.R
index 02ff0a8..944daf1 100644
--- a/R/domain.R
+++ b/R/domain.R
@@ -16,19 +16,11 @@
                       )
 )
 
-#' @name Domain
-#'
-#' @exportClass Domain
 setOldClass(c("Domain", "R6"))
 
-#' Mesh 
-#'
-#' @name Mesh
-#'
-#' @exportClass Mesh
 setOldClass(c("Mesh","Domain"))
 
-#' S4 class representing a spatial (spatio-temporal) domain
+#' R6 class representing a spatial (spatio-temporal) domain
 #' 
 #' @param x could be a \code{pslg} returned by \code{\link[RTriangle]{triangulate}}, 
 #' a \code{sfc} returned by \code{\link[sf]{st_geometry}} or a named list containing:
@@ -294,7 +286,7 @@ setMethod("build_mesh", signature=c("Domain", "numeric", "numeric"),
 #'
 #' @param op1 A mesh object created by \code{Mesh}.
 #' @param op2 A numeric vector.
-#' @return An S4 object representing a spatio-temporal domain.
+#' @return A R6 object representing a spatio-temporal domain.
 #' @rdname Domain_times_vector
 #' @export 
 setGeneric("%X%", function(op1, op2) standardGeneric("%X%"))
diff --git a/R/function_space.R b/R/function_space.R
index b6401e5..41238d6 100644
--- a/R/function_space.R
+++ b/R/function_space.R
@@ -32,11 +32,6 @@
   )
 )
 
-#' Finite Element Basis 
-#'
-#' @name Basis
-#'
-#' @exportClass BasisFunction
 setOldClass(c("BasisFunction", "R6"))
 
 setGeneric("BasisFunction", function(mesh, fe_order) standardGeneric("BasisFunction"))
@@ -74,18 +69,13 @@ setMethod("BasisFunction", signature = signature("Mesh","integer"),
                               )
 )
 
-#' Finite Element Functional Space
-#'
-#' @name FunctionSpace
-#'
-#' @exportClass FunctionSpace
 setOldClass(c("FunctionSpace", "R6"))
 
 #' Create FunctionSpace object
 #'
 #' @param mesh A mesh object created by \code{Mesh}:
 #' @param fe_order Either '1' or '2'. It specifies the finite element order.
-#' @return An S4 object representing a Function Space.
+#' @return A R6 object representing a Function Space.
 #' @export 
 #' @rdname FunctionSpace
 setGeneric("FunctionSpace", function(mesh, fe_order) standardGeneric("FunctionSpace"))
@@ -169,11 +159,6 @@ setMethod("FunctionSpace", signature = c(mesh="Mesh", fe_order="missing"),
                         )
 )
 
-#' Finite Element Function
-#'
-#' @name Function
-#'
-#' @exportClass Function
 setOldClass(c("Function", "R6"))
 
 ## constructor
@@ -181,7 +166,7 @@ setOldClass(c("Function", "R6"))
 #' Create Function object
 #'
 #' @param FunctionSpace a FunctionSpace object created by \code{FunctionSpace}:
-#' @return An S4 object representing a Function belonging to the FunctionSpace passed as parameter.
+#' @return A R6 object representing a Function belonging to the FunctionSpace passed as parameter.
 #' @export 
 #' @examples
 #' \dontrun{
diff --git a/R/modules.R b/R/modules.R
index 6d7b8d7..fc405fc 100644
--- a/R/modules.R
+++ b/R/modules.R
@@ -4,6 +4,7 @@
 #' @importFrom Rcpp evalCpp
 #' @import plotly
 #' @import Matrix
+#' @import R6
 NULL
 
 ## load required modules
diff --git a/R/operators.R b/R/operators.R
index 522a043..1d6812e 100644
--- a/R/operators.R
+++ b/R/operators.R
@@ -15,11 +15,6 @@
   )
 )
 
-#' Function gradient
-#'
-#' @name grad
-#'
-#' @exportClass FunctionGrad
 setOldClass(c("FunctionGrad", "R6"))
 
 ## take gradient of Function
@@ -27,7 +22,7 @@ setOldClass(c("FunctionGrad", "R6"))
 #' compute gradient of Function
 #'
 #' @param Function a Function created by \code{Function}:
-#' @return An S4 object representing the gradient of the Function passed as parameter.
+#' @return A R6 object representing the gradient of the Function passed as parameter.
 #' @export 
 #' @rdname grad
 #' @examples
@@ -46,7 +41,7 @@ setMethod("grad", signature(Function = "Function"), function(Function) {
   .FunctionGradCtr$new(Function = Function)
 })
 
-# setMethod("*", signature=c(e1="matrix", e2="FunctionGrad")
+# #`*.FunctionGrad` <- function(e1, e2){ setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"),
 
 #' product overload for FunctionGradObejct
 #'
@@ -60,7 +55,7 @@ setMethod("grad", signature(Function = "Function"), function(Function) {
 }
 
 ## base class for differential operators
-.DiffOpCtr <- R6Class("DiffOpObject",
+.DifferentialOpCtr <- R6Class("DifferentialOp",
     private = list(
       ## Function to which the operator is applied
       Function_ = "Function"
@@ -80,25 +75,20 @@ setMethod("grad", signature(Function = "Function"), function(Function) {
     )
 )
 
-#' R6 class representing a differential operator
-#' 
-#' @name Differential Operator
-#'
-#' @exportClass DiffOpObject
-setOldClass(c("DiffOpObject", "R6"))
+setOldClass(c("DifferentialOp", "R6"))
 
-##setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
+##setMethod("+", signature = c(e1="DifferentialOp", e2="DifferentialOp"),
+## `+.DifferentialOp` <- function(e1, e2){
 
 ## sum of differential operators
 
-#' plus operator overload for DiffOpObject
+#' Sum of DifferentialOp
 #'
-#' @param e1 a DiffOpObject.
-#' @param e2 a DiffOpObject.
-#' @return A S4 object representing the sum of two differential operators.
-#' @name DifferentialOps
+#' @param e1 a \code{DifferentialOp}.
+#' @param e2 a \code{DifferentialOp}.
+#' @return A R6 class representing the sum of two differential operators.
 #' @export 
-`+.DiffOpObject` <- function(e1, e2){
+`+.DifferentialOp` <- function(e1, e2){
             if (!identical(e1$Function(), e2$Function())) {
               stop("operator arguments must be the same")
             }
@@ -106,30 +96,32 @@ setOldClass(c("DiffOpObject", "R6"))
             if (any(duplicated(c(e1$tokens, e2$tokens)))) {
               stop("duplicated operator")
             }
-            .DiffOpCtr$new(
+            .DifferentialOpCtr$new(
               tokens = c(e1$tokens, e2$tokens),
               params = c(e1$params, e2$params),
               Function = e1$Function()
             )
 }
 
-# setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
+# setMethod("-", signature = c(e1="DifferentialOp", e2="DifferentialOp"),
+# `-.DifferentialOp` <- function(e1, e2){
 
-#' difference of differential operators
+#' Difference of differential operators 
 #'
-#' @param e1 a DiffOpObject.
-#' @param e2 a DiffOpObject.
-#' @return A S4 object representing the difference of two differential operators.
-#' @rdname minus_DiffOb_op
-#' @name DifferentialOps
+#' @param e1 a \code{DifferentialOp}.
+#' @param e2 a \code{DifferentialOp} or missing.
+#' @return A R6 class representing the difference of two differential operators.
 #' @export 
-`-.DiffOpObject` <- function(e1, e2){
+`-.DifferentialOp` <-function(e1 ,e2){
             if(missing(e2)){
               e1$params[[1]] <- -e1$params[[1]]
               return(e1)
             }else if (!identical(e1$Function(), e2$Function())) {
               stop("operator arguments must be the same")
             }
+            if (!identical(e1$Function(), e2$Function())) 
+              stop("operator arguments must be the same")
+            
             ## check for duplicated operator tokens
             if (any(duplicated(c(e1$tokens, e2$tokens)))) {
               stop("duplicated operator")
@@ -137,37 +129,44 @@ setOldClass(c("DiffOpObject", "R6"))
             for(i in 1:length(e2$params)){
               e2$params[[i]] = -e2$params[[i]] 
             }
-            .DiffOpCtr$new(
+            .DifferentialOpCtr$new(
               tokens = c(e1$tokens, e2$tokens),
               params = c(e1$params, e2$params),
               Function = e1$Function()
             )
-  }
+}
 
-#setMethod("-", signature =c(e1 = "DiffOpObject", e2 = "missing"),
+#setMethod("-", signature =c(e1 = "DifferentialOp", e2 = "missing"),
 
 # minus (unary) operator for differential operators
 # 
-# @param e1 a DiffOpObject.
-# @return A S4 object differential operator whose paramter has changed sign.
+# @param e1 a DifferentialOp.
+# @return A R6 class differential operator whose paramter has changed sign.
+# @name DifferentialOps
 # @export
+# setMethod("-", signature =c(e1 = "DifferentialOp", e2 = "missing"),
+#           function(e1 ,e2){
+#             e1$params[[1]] <- -e1$params[[1]]
+#             e1
+# })
 
 # @rdname minus_DiffOb_op
-# `-.DiffOpObject` <- function(e1){
+# `-.DifferentialOp` <- function(e1){
 #             e1$params[[1]] <- -e1$params[[1]]
 #             e1
 # }
 
-# setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"),
+# setMethod("*", signature=c(e1="numeric", e2="DifferentialOp"),
+# `*.DifferentialOp` <- function(e1,e2){  
 
 #' product by scalar for differential operators
 #'
 #' @param e1 a numeric.
-#' @param e2 a DiffOpObject.
-#' @return A S4 object representing the product by scalar for a differential operator.
-#' @name DifferentialOps
+#' @param e2 a DifferentialOp.
+#' @return A R6 class representing the product by scalar for a differential operator.
+#' @name times_DifferentialOps
 #' @export 
-`*.DiffOpObject` <- function(e1,e2){  
+`*.DifferentialOp` <- function(e1,e2){
             if (!is.numeric(e1)) stop("bad product")
             e2$params[[1]] <- e1*e2$params[[1]]
             e2
@@ -175,15 +174,10 @@ setOldClass(c("DiffOpObject", "R6"))
   
 ## diffusion term
 .DiffusionCtr <- R6Class("DiffusionOperator",
-    inherit = .DiffOpCtr
+    inherit = .DifferentialOpCtr
 )
 
-#' Diffusion Operator
-#'
-#' @name DifferentialOperator
-#'
-#' @exportClass DiffusionOperator
-setOldClass(c("DiffusionOperator", "DiffOpObject"))
+setOldClass(c("DiffusionOperator", "DifferentialOp"))
 
 ## laplace() returns a special operator for the case of
 ## isotropic  and stationary diffusion
@@ -191,7 +185,7 @@ setOldClass(c("DiffusionOperator", "DiffOpObject"))
 #' laplace operator for Function class
 #'
 #' @param Function a Function.
-#' @return A S4 object representing the laplace operator applied to the function passed as parameter.
+#' @return A R6 class representing the laplace operator applied to the function passed as parameter.
 #' @rdname DiffusionOperator
 #' @export 
 #' @examples
@@ -219,7 +213,7 @@ laplace <- function(Function) {
 #' divergence operator FunctionGrad
 #'
 #' @param Function a Function.
-#' @return A S4 object representing the diffusion term of a second order linear differential operator.
+#' @return A R6 class representing the diffusion term of a second order linear differential operator.
 #' @rdname DiffusionOperator
 #' @export
 #' @examples
@@ -247,21 +241,16 @@ div <- function(Function) {
 
 ## transport term
 .TransportCtr <- R6Class("TransportOperator",
-    inherit = .DiffOpCtr
+    inherit = .DifferentialOpCtr
 )
 
-#' Transport Operator
-#'
-#' @name TransportOperator
-#'
-#' @exportClass TransportOperator
-setOldClass(c("TransportOperator", "DiffOpObject"))
+setOldClass(c("TransportOperator", "DifferentialOp"))
 
 #' dot product between vector and FunctionGrad
 #'
 #' @param op1 a numeric vector.
 #' @param op2 a FunctionGrad.
-#' @return A S4 object representing the advection term of a second order linear differential operator.
+#' @return A R6 class representing the advection term of a second order linear differential operator.
 #' @rdname TransportOperator
 #' @export
 #' @examples
@@ -288,23 +277,19 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"),
 
 ## reaction term
 .ReactionCtr <- R6Class("ReactionOperator",
-    inherit = .DiffOpCtr
+    inherit = .DifferentialOpCtr
 )
 
-#' Reaction Operator
-#'
-#' @name ReactionOperator
-#'
-#' @exportClass ReactionOperator
-setOldClass(c("ReactionOperator", "DiffOpObject"))
+setOldClass(c("ReactionOperator", "DifferentialOp"))
 
 # setMethod("*", signature = c(e1="numeric", e2="Function"),
-          
+# `*.Function` <- function(e1,e2){
+
 #' product by scalar for FunctionObejct
 #'
 #' @param e1 a numeric.
 #' @param e2 a FunctioObject created by \code{Function}.
-#' @return A S4 object representing the reaction term of a second order linear differential operator.
+#' @return A R6 class representing the reaction term of a second order linear differential operator.
 #' @rdname ReactionOperator
 #' @export 
 #' @examples
@@ -326,13 +311,10 @@ setOldClass(c("ReactionOperator", "DiffOpObject"))
 }
 
 .TimeDerivativeCtr <- R6Class("TimeDerivative",
-  inherit = .DiffOpCtr
+  inherit = .DifferentialOpCtr
 )
 
-#' @name TimeDerivative
-#'
-#' @exportClass TimeDerivative
-setOldClass(c("TimeDerivative", "DiffOpObject"))
+setOldClass(c("TimeDerivative", "DifferentialOp"))
 
 # overloading stats::dt function
 
@@ -341,7 +323,7 @@ setOldClass(c("TimeDerivative", "DiffOpObject"))
 #' @param x a Function.
 #' @param df missing.
 #' @param ncp missing.
-#' @return A S4 object representing the time derivative of a Function.
+#' @return A R6 class representing the time derivative of a Function.
 #' @rdname TimeDerivative
 #' @export
 #' @examples
diff --git a/R/pde.R b/R/pde.R
index c6e1144..6d402c7 100644
--- a/R/pde.R
+++ b/R/pde.R
@@ -8,7 +8,7 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3)
     is_initialCondition_set = "logical",
     is_parabolic ="logical",          
     cpp_handler_ = "ANY",
-    DifferentialOp_ = "DiffOpObject"
+    DifferentialOp_ = "DifferentialOp"
   ),
   public = list(
     initialize = function(cpp_handler, DifferentialOp, is_parabolic, is_dirichletBC_set, is_initialCondition_set){
@@ -40,26 +40,26 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3)
     dofs_coordinates = function(){
       private$cpp_handler_$dofs_coordinates()
     },
-    set_boundary_condition = function(fun, type="dirichlet", on=NULL){
+    set_boundary_condition = function(boundary_condition, type="dirichlet", on=NULL){
       if(!any(type == c("dirichlet", "Dirichlet", "d"))) 
         stop("Only Dirichlet boundary condtions allowed.")
       dirichletBC_ <- NULL
-      if(any(typeof(fun) == c("function", "closure"))){
+      if(any(typeof(boundary_condition) == c("function", "closure"))){
         if(!private$is_parabolic){
-          dirichletBC_ <- as.matrix(fun(private$cpp_handler_$dofs_coordinates()))
+          dirichletBC_ <- as.matrix(boundary_condition(private$cpp_handler_$dofs_coordinates()))
         }else{
           times <- private$DifferentialOp_$Function()$FunctionSpace()$mesh()$time_nodes()
-          dirichletBC_ <- fun(private$cpp_handler_$dofs_coordinates(), times)
+          dirichletBC_ <- boundary_condition(private$cpp_handler_$dofs_coordinates(), times)
         }
-      }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){
-        if(nrow(as.matrix(fun)) == nrow(private$cpp_handler_$dofs_coordinates())){
-          dirichletBC_ <- fun  
-        }else if (nrow(as.matrix(fun)) == 1L){
+      }else if(any(typeof(boundary_condition) == c("matrix","vector", "numeric" ,"double"))){
+        if(nrow(as.matrix(boundary_condition)) == nrow(private$cpp_handler_$dofs_coordinates())){
+          dirichletBC_ <- boundary_condition  
+        }else if (nrow(as.matrix(boundary_condition)) == 1L){
           if(!private$is_parabolic)
-            dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler_$dofs_coordinates()), ncol=1)
+            dirichletBC_ <- matrix(boundary_condition, nrow=nrow(private$cpp_handler_$dofs_coordinates()), ncol=1)
           else{
             times <- private$DifferentialOp_$Function()$FunctionSpace()$mesh()$time_nodes()
-            dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler_$dofs_coordinates()), 
+            dirichletBC_ <- matrix(boundary_condition, nrow=nrow(private$cpp_handler_$dofs_coordinates()), 
                                         ncol=length(times))
           }
         }   
@@ -68,17 +68,17 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3)
       private$cpp_handler_$set_dirichlet_bc(dirichletBC_)
       invisible(self)
     },
-    set_initial_condition = function(fun){
+    set_initial_condition = function(initial_condition){
       if(!private$is_parabolic)
         stop("Cannot set initial condition for elliptic problem.")
       private$is_initialCondition_set <- TRUE
-      if(typeof(fun) == "closure"){ 
-        private$cpp_handler_$set_initial_condition(fun(private$cpp_handler_$dofs_coordinates()))
-      }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){
-        if(nrow(as.matrix(fun)) == nrow(private$cpp_handler_$dofs_coordinates())){
-          private$cpp_handler_$set_initial_condition(fun)
-        }else if(nrow(as.matrix(fun)) == 1L){
-          private$cpp_handler_$set_initial_condition(matrix(fun, nrow=nrow(private$cpp_handler_$dofs_coordinates()),
+      if(typeof(initial_condition) == "closure"){ 
+        private$cpp_handler_$set_initial_condition(initial_condition(private$cpp_handler_$dofs_coordinates()))
+      }else if(any(typeof(initial_condition) == c("matrix","vector", "numeric" ,"double"))){
+        if(nrow(as.matrix(initial_condition)) == nrow(private$cpp_handler_$dofs_coordinates())){
+          private$cpp_handler_$set_initial_condition(initial_condition)
+        }else if(nrow(as.matrix(initial_condition)) == 1L){
+          private$cpp_handler_$set_initial_condition(matrix(initial_condition, nrow=nrow(private$cpp_handler_$dofs_coordinates()),
                                                  ncol=1))   
         }
       }
@@ -93,9 +93,6 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3)
   )
 )
 
-#' @name pde
-#'
-#' @exportClass Pde
 setOldClass(c("Pde", "R6"))
 
 ## infers the type of a pde
@@ -154,15 +151,15 @@ make_pde <- function(DifferentialOp) {
 #' A PDEs object
 #'
 #' @param DifferentialOp a differential operator.
-#' @param f a standard R function representing the forcing term of the PDE or a numeric value, in case of constant forcing term.
+#' @param forcing a standard R function representing the forcing term of the PDE or a numeric value, in case of constant forcing term.
 #' @return A R6 class representing a PDE.
 #' @rdname pde
 #' @export 
-setGeneric("Pde", function(DifferentialOp, f) standardGeneric("Pde"))
+setGeneric("Pde", function(DifferentialOp, forcing) standardGeneric("Pde"))
 
 #' @rdname pde
-setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"),
-          function(DifferentialOp,f){
+setMethod("Pde", signature=c(DifferentialOp="DifferentialOp", forcing="function"),
+          function(DifferentialOp,forcing){
             
             pde_type = extract_pde_type(DifferentialOp)
             cpp_handler <- make_pde(DifferentialOp)
@@ -171,9 +168,9 @@ setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"),
             ## evaluate forcing term on quadrature nodes
             if(pde_type == pde_type_list$parabolic) {
               times <- DifferentialOp$Function()$FunctionSpace()$mesh()$time_nodes()
-              cpp_handler$set_forcing(as.matrix(f(quad_nodes, times)))
+              cpp_handler$set_forcing(as.matrix(forcing(quad_nodes, times)))
             }else{
-              cpp_handler$set_forcing(as.matrix(f(quad_nodes)))
+              cpp_handler$set_forcing(as.matrix(forcing(quad_nodes)))
             }
             
             ## initialize solver 
@@ -190,20 +187,20 @@ setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"),
 })
 
 #' @rdname pde
-setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="numeric"),
-          function(DifferentialOp,f){
+setMethod("Pde", signature=c(DifferentialOp="DifferentialOp", forcing="numeric"),
+          function(DifferentialOp,forcing){
             pde_type = extract_pde_type(DifferentialOp)
             
-            fun <- NULL
+            forcing_ <- NULL
             if(pde_type!=pde_type_list$parabolic){
-              fun <- function(points){
-                return( matrix(f, nrow=nrow(points), ncol=1))
+              forcing_ <- function(points){
+                return( matrix(forcing, nrow=nrow(points), ncol=1))
               }
             }else{
-              fun <- function(points, times){
-                return( matrix(f, nrow=nrow(points), ncol=length(times)))
+              forcing_ <- function(points, times){
+                return( matrix(forcing, nrow=nrow(points), ncol=length(times)))
               }  
             }
             
-            return(Pde(DifferentialOp,fun))
+            return(Pde(DifferentialOp,forcing_))
 })
diff --git a/R/sf.R b/R/sf.R
index 24db96a..a7e4f4c 100644
--- a/R/sf.R
+++ b/R/sf.R
@@ -1,13 +1,12 @@
-#' sf methods for Domain and Mesh
+#' Integration with sf package
 #'
-#' \code{\link[sf]{sf}} methods for \code{Domain} and \code{Mesh} objects.
+#' \code{\link[sf]{sf}} methods for \code{Domain} and \code{Mesh} classes.
 ##
 #' @param x An object of class \code{Domain} or \code{Mesh}.
-#' @param ... Arguments passed on the corresponding \code{sf} function.
+#' @param ... Arguments passed to the corresponding \code{sf} method
 #' @param value The value to be assigned. See the documentation of the
 #' corresponding sf function for details.
-#' @return The \code{femR} method for \code{\link[sf]{st_as_sf}} returns
-#' ............ \code{\link[sf]{sf}}.
+#' @return All methods return the same type of object of the corresponding \code{\link[sf]{sf}} method.
 #' @importFrom sf st_as_sf
 #' @importFrom sf st_polygon
 #' @importFrom sf st_as_sfc
diff --git a/R/utils.R b/R/utils.R
index 9a2968f..5aaf514 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -10,8 +10,3 @@ setMethod("set_private", signature = c("R6", "character"),
           function(x, attribute, value){
             invisible(x$.__enclos_env__$private[[attribute]] <- value)
 })
-
-# extract_private <- function(x){
-#   if(!is(x,"R6")) stop("Input parameter must be an R6 class.")
-#   x$.__enclos_env__$private
-# }
diff --git a/man/Basis.Rd b/man/Basis.Rd
deleted file mode 100644
index ddd13f4..0000000
--- a/man/Basis.Rd
+++ /dev/null
@@ -1,8 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/function_space.R
-\name{Basis}
-\alias{Basis}
-\title{Finite Element Basis}
-\description{
-Finite Element Basis
-}
diff --git a/man/Differential-Operator.Rd b/man/Differential-Operator.Rd
deleted file mode 100644
index 31149b3..0000000
--- a/man/Differential-Operator.Rd
+++ /dev/null
@@ -1,8 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{Differential Operator}
-\alias{Differential Operator}
-\title{R6 class representing a differential operator}
-\description{
-R6 class representing a differential operator
-}
diff --git a/man/DifferentialOperator.Rd b/man/DifferentialOperator.Rd
deleted file mode 100644
index 18948ec..0000000
--- a/man/DifferentialOperator.Rd
+++ /dev/null
@@ -1,8 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{DifferentialOperator}
-\alias{DifferentialOperator}
-\title{Diffusion Operator}
-\description{
-Diffusion Operator
-}
diff --git a/man/DifferentialOps.Rd b/man/DifferentialOps.Rd
deleted file mode 100644
index 7b5962e..0000000
--- a/man/DifferentialOps.Rd
+++ /dev/null
@@ -1,27 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{DifferentialOps}
-\alias{DifferentialOps}
-\alias{+.DiffOpObject}
-\alias{*.DiffOpObject}
-\title{plus operator overload for DiffOpObject}
-\usage{
-\method{+}{DiffOpObject}(e1, e2)
-
-\method{*}{DiffOpObject}(e1, e2)
-}
-\arguments{
-\item{e1}{a numeric.}
-
-\item{e2}{a DiffOpObject.}
-}
-\value{
-A S4 object representing the sum of two differential operators.
-
-A S4 object representing the product by scalar for a differential operator.
-}
-\description{
-plus operator overload for DiffOpObject
-
-product by scalar for differential operators
-}
diff --git a/man/minus_DiffOb_op.Rd b/man/minus_DiffOb_op.Rd
deleted file mode 100644
index dc352af..0000000
--- a/man/minus_DiffOb_op.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{DifferentialOps}
-\alias{DifferentialOps}
-\alias{-.DiffOpObject}
-\title{difference of differential operators}
-\usage{
-\method{-}{DiffOpObject}(e1, e2)
-}
-\arguments{
-\item{e1}{a DiffOpObject.}
-
-\item{e2}{a DiffOpObject.}
-}
-\value{
-A S4 object representing the difference of two differential operators.
-}
-\description{
-difference of differential operators
-}
diff --git a/src/core b/src/core
index 8fa5175..182556d 160000
--- a/src/core
+++ b/src/core
@@ -1 +1 @@
-Subproject commit 8fa5175f22f6479ec62723832906ae2e71b74019
+Subproject commit 182556d3e1c93793f88d31b15f9eed6f6d07967d
diff --git a/tests/parabolic.R b/tests/parabolic.R
index c9065cb..449530b 100644
--- a/tests/parabolic.R
+++ b/tests/parabolic.R
@@ -54,7 +54,7 @@ u0 <- function(points){
 }
   
 ## create pde
-pde <- Pde(L, f)
+pde <- Pde(L, 0.)
 pde$set_boundary_condition(g)
 pde$set_initial_condition(u0)
 
diff --git a/tests/pde.2.R b/tests/pde.2.R
index 5af2f04..7f177a8 100644
--- a/tests/pde.2.R
+++ b/tests/pde.2.R
@@ -28,7 +28,7 @@ g <- function(points){
 }
 ## Pde constructor
 pde <- Pde(L,u)
-pde$set_boundary_condition(fun=0.,
+pde$set_boundary_condition(boundary_condition =0.,
                            type = "dirichlet")
 
 ## solve problem
diff --git a/vignettes/Introduction.Rmd b/vignettes/Introduction.Rmd
index 102c9b4..cd62d57 100644
--- a/vignettes/Introduction.Rmd
+++ b/vignettes/Introduction.Rmd
@@ -106,7 +106,7 @@ g <- function(points){
 pde <- Pde(Lu, f)
 # setting boundary conditions
 pde$set_boundary_condition(type="dirichlet",
-                           fun=g)
+                           boundary_condition=g)
 ## 7. computing the discrete solution
 pde$solve()
 
@@ -255,7 +255,7 @@ pde <- Pde(Lu, f)
 
 ## set Dirichlet boundary conditions
 pde$set_boundary_condition(type="dirichlet",
-                           fun=g)
+                           boundary_condition=g)
 ```
 We can compute the discrete solution of the Poisson problem calling the `solve` method:
 ```{r Solving the problem}
@@ -342,7 +342,7 @@ u0 <- function(points){
 pde <- Pde(Lu, f)
 # set initial condition and dirichlet BC
 pde$set_initial_condition(u0)
-pde$set_boundary_condition(fun= 0.0, type="dirichlet")
+pde$set_boundary_condition(boundary_condition= 0.0, type="dirichlet")
 ## solve problem
 pde$solve()
 ## plot solution
diff --git a/vignettes/LakeComo.Rmd b/vignettes/LakeComo.Rmd
index 2bac927..3f9cfab 100644
--- a/vignettes/LakeComo.Rmd
+++ b/vignettes/LakeComo.Rmd
@@ -87,15 +87,10 @@ The upcoming sections will illustrate how to obtain a discrete solution of the p
 library(sf, quietly = T)
 library(mapview, quietly = T)
 library(femR, quietly = T)
-library(rmapshaper, quietly = T)
 
-lake_bd <- read_sf("data/deims_sites_boundariesPolygon.shp")
-lake_como_boundary <- rmapshaper::ms_simplify(lake_bd, keep= 0.5, 
-                                        keep_shapes=T)
+data("lake_como_boundary", package="femR")
 st_crs(lake_como_boundary) <- 4326
-
 ```
-
 The `femR` package includes utilities  to read a `sfc` (Simple Features Collection) object, such as the boundary of the Lake, and generating meshes by specifying the `maximum_area` of each triangle and the `minimum_angle` between adjacent sides of the triangles.
 ```{r}
 # create the physical domain
@@ -149,7 +144,7 @@ Finally, we build a `Pde` object passing the differential operator `L`, as first
 # build PDE object
 pde <- Pde(Lu, 0.)
 # set boundary conditions
-pde$set_boundary_condition(fun= g, type= "dirichlet")
+pde$set_boundary_condition(boundary_condition= g, type= "dirichlet")
 ```
 We can compute the discrete solution of the problem calling the `solve` method:
 ```{r}
@@ -264,7 +259,7 @@ pde <- Pde(Lu, 0.)
 # set initial conditions
 pde$set_initial_condition(u0)
 # set boundary conditions
-pde$set_boundary_condition(fun= g, type= "dirichlet")
+pde$set_boundary_condition(boundary_condition= g, type= "dirichlet")
 # compute the discrete solution
 pde$solve()
 ```