diff --git a/DESCRIPTION b/DESCRIPTION
index cd786a0..c6e8aed 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -15,7 +15,7 @@ Authors@R: c( person("Aldo", "Clemente", role = c("aut", "cre"),
 	                                        email = "laura.sangalli@polimi.it"))
 Maintainer: Aldo Clemente <aldo.clemente@polimi.it>
 Depends: 
-    R (>= 3.5.0),
+    R (>= 3.5.0), R6,
     plotly, Matrix, methods,
     RTriangle, sf
 Description: Solving partial differential equations relying on the Finite Element Method.
diff --git a/NAMESPACE b/NAMESPACE
index 106ec36..4299a24 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method("*",DiffOpObject)
+S3method("*",Function)
+S3method("*",FunctionGrad)
+S3method("+",DiffOpObject)
+S3method("-",DiffOpObject)
 S3method("st_crs<-",Domain)
 S3method(plot,Function)
 S3method(plot,Mesh)
@@ -8,6 +13,8 @@ S3method(st_as_sfc,Domain)
 S3method(st_as_sfc,Mesh)
 S3method(st_crs,Domain)
 export("%X%")
+export(.FunctionSpaceCtr)
+export(.MeshCtr)
 export(Domain)
 export(Function)
 export(FunctionSpace)
@@ -19,9 +26,18 @@ export(dot)
 export(grad)
 export(laplace)
 export(read_mesh)
-exportMethods("*")
-exportMethods("+")
-exportMethods("-")
+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)
diff --git a/R/domain.R b/R/domain.R
index a16e768..c69fcc4 100644
--- a/R/domain.R
+++ b/R/domain.R
@@ -1,12 +1,29 @@
-.DomainCtr <- setRefClass("Domain", 
-                          fields=c(
-                            geometry = "ANY",
-                            coords = "data.frame",         # (x,y,boundary)
-                            time_interval = "numeric",
-                            crs ="ANY"
-                           )
+.DomainCtr <- R6Class("Domain", 
+                      public = list(
+                        geometry = "ANY",
+                        coords = "data.frame",         # (x,y,boundary)
+                        time_interval = vector(mode="numeric", length = 0),
+                        crs = "ANY",
+                        initialize = function(geometry, coords, crs){
+                                      self$geometry <- geometry 
+                                      self$coords <- coords 
+                                      self$crs <- crs
+                                    }
+                      )
 )
 
+#' @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
 #' 
 #' @param x could be a \code{pslg} returned by \code{\link[RTriangle]{triangulate}}, 
@@ -36,10 +53,11 @@ setMethod("Domain", signature = "list", function(x){
     
     coords <- data.frame(x=x$nodes[,1], y=x$nodes[,2], boundary=x$nodes_boundary)
   crs = NA
-  .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
-             coords=coords, crs = crs)
+  .DomainCtr$new(geometry = geometry, coords=coords, crs = crs)
 })
 
+setOldClass("pslg")
+
 #' @rdname Domain
 setMethod("Domain", signature = "pslg", function(x){
   nodes <- x$P
@@ -62,8 +80,7 @@ setMethod("Domain", signature = "pslg", function(x){
   
   coords <- data.frame(x=nodes[,1], y=nodes[,2], boundary=nodes_boundary)
   
-  .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0),
-             coords=coords, crs = NA)
+  .DomainCtr$new(geometry = geometry, coords=coords, crs = NA)
 })
 
 #' @rdname Domain
@@ -206,8 +223,7 @@ setMethod("Domain", signature="sfc", function(x){
   
   crs <- NA
   if(!is.na(st_crs(x))) crs <- st_crs(x)$input
-      .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
-                coords=coords, crs = crs)
+      .DomainCtr$new(geometry = geometry, coords=coords, crs = crs)
 })
 
 # setMethod("Domain", signature = "sfc", function(x){
@@ -222,7 +238,7 @@ setMethod("Domain", signature="sfc", function(x){
 #                    nodes_group = nodes_group, edges_group = edges_group)
 #   crs <- NA_crs_
 #   if(!is.na(st_crs(x))) crs <- st_crs(x)
-#   .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
+#   .DomainCtr$new(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
 #              coords=geometry$nodes, crs = crs)
 # })
 
diff --git a/R/function_space.R b/R/function_space.R
index d5b228c..9f7174f 100644
--- a/R/function_space.R
+++ b/R/function_space.R
@@ -1,49 +1,68 @@
-
-.BasisFunctionCtr <- setRefClass(
-  Class = "BasisFunction",
-  fields = c(
-    cpp_handler = "ANY"  ## pointer to cpp backend
-  ),
-  methods = c(
-    ## evaluates the basis system on the given set of locations
-    eval = function(locations) {
-      return(cpp_handler$eval(0L, locations))
-    },
-    ## integrates a function expressed as basis expansion with respect to this basis system over the
-    ## whole domain of definition
-    integrate = function(func) {
-      if (!(is.function(func) || is.vector(func))) stop("invalid argument, should be either a function or a vector")
-      ## required dof coordinates...
-      c <- func
-      if (is.function(func)) { ## recover fem expansion
-        c <- func(mesh$cpp_handler$nodes())
-      }
-      return(cpp_handler$integrate(c))
-    }
+.BasisFunctionCtr <- R6Class("BasisFunction",
+                              public = list(
+                                cpp_handler = "ANY",  ## pointer to cpp backend
+                                mesh = "Mesh",
+                                initialize = function(cpp_handler, mesh){
+                                  self$cpp_handler = cpp_handler
+                                  self$mesh = mesh
+                                },
+                                ## evaluates the basis system on the given set of locations
+                                eval = function(locations) {
+                                  return(self$cpp_handler$eval(0L, locations))
+                                },
+                                ## integrates a function expressed as basis expansion with respect to this basis system over the
+                                ## whole domain of definition
+                                integrate = function(func) {
+                                  if (!(is.function(func) || is.vector(func))) stop("invalid argument, should be either a function or a vector")
+                                  ## required dof coordinates...
+                                  c <- func
+                                  if (is.function(func)) { ## recover fem expansion
+                                    c <- func(self$mesh$cpp_handler$nodes())
+                                  }
+                                  return(cpp_handler$integrate(c))
+                                }
   )
 )
 
+#' Finite Element Basis 
+#'
+#' @name Basis
+#'
+#' @exportClass BasisFunction
+setOldClass(c("BasisFunction", "R6"))
+
 setGeneric("BasisFunction", function(mesh, fe_order) standardGeneric("BasisFunction"))
 
 setMethod("BasisFunction", signature = signature("Mesh","integer"),
           function(mesh, fe_order){
-            basis_function <- .BasisFunctionCtr(cpp_handler = 
+            basis_function <- .BasisFunctionCtr$new(cpp_handler = 
                                                   new(eval(parse(text = paste("cpp_lagrange_basis_2d_fe", 
-                                                                              as.character(fe_order), sep = ""))), mesh$cpp_handler, 0))
+                                                                              as.character(fe_order), sep = ""))), mesh$cpp_handler, 0),
+                                                  mesh=mesh)
 })
 
-.FunctionSpaceCtr <- setRefClass(
-  Class = "FunctionSpaceObject",
-  fields = c(
-    mesh = "ANY",
-    basis_function = "BasisFunction",     ## cpp backend
-    fe_order = "integer"                  ## this is specific for fem, must be generalized
-  ),
-  methods = c(
-    get_basis = function() { return(basis_function) }
-  )
+#' @export 
+.FunctionSpaceCtr <- R6Class("FunctionSpace",
+                              public = list(
+                                mesh = "ANY",
+                                basis_function = "BasisFunction",     ## cpp backend
+                                fe_order = "integer",                  ## this is specific for fem, must be generalized
+                                initialize = function(mesh, basis_function, fe_order){
+                                  self$mesh <- mesh 
+                                  self$basis_function <- basis_function
+                                  self$fe_order <- fe_order
+                                },
+                                get_basis = function() { return(self$basis_function) }
+                              )
 )
 
+#' Finite Element Functional Space
+#'
+#' @name FunctionSpace
+#'
+#' @exportClass FunctionSpace
+setOldClass(c("FunctionSpace", "R6"))
+
 #' Create FunctionSpace object
 #'
 #' @param mesh A mesh object created by \code{Mesh}:
@@ -51,7 +70,7 @@ setMethod("BasisFunction", signature = signature("Mesh","integer"),
 #' @return An S4 object representing a Function Space.
 #' @export 
 #' @rdname FunctionSpace
-setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpace"))
+setGeneric("FunctionSpace", function(mesh, fe_order) standardGeneric("FunctionSpace"))
 
 #' @rdname FunctionSpace
 #' @examples
@@ -62,11 +81,11 @@ setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpa
 #' Vh <- FunctionSpace(mesh = mesh, fe_order = 1)
 #' }
 setMethod("FunctionSpace",
-          signature = c(mesh = "ANY", fe_order = "numeric"),
+          signature = c(mesh = "Mesh", fe_order = "numeric"),
           function(mesh, fe_order) {
-            .FunctionSpaceCtr(
-              basis_function = BasisFunction(mesh,as.integer(fe_order)),
+            .FunctionSpaceCtr$new(
               mesh = mesh,
+              basis_function = BasisFunction(mesh, as.integer(fe_order)),
               fe_order = as.integer(fe_order)
             )
           }
@@ -80,46 +99,55 @@ setMethod("FunctionSpace",
 #' mesh <- Mesh(unit_square)
 #' Vh <- FunctionSpace(mesh)
 #' }
-setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
+setMethod("FunctionSpace", signature = c(mesh="Mesh", fe_order="missing"),
           function(mesh){
-              return(.FunctionSpaceCtr(mesh=mesh, fe_order=1L))
+              return(.FunctionSpaceCtr$new(mesh=mesh, 
+                                           basis_function = BasisFunction(mesh, 1L), 
+                                           fe_order=1L))
 })
 
 ## finite element function
-.FunctionCtr <- setRefClass(
-  Class = "Function",
-  fields = c(
-    FunctionSpace  = "ANY", 
-    coeff = "matrix",
-    pde = "ANY"
-  ),
-  methods = list(
-    eval = function(X) {
-      M = dim(FunctionSpace$mesh$get_nodes())[2]
-      if(is.vector(X)) X <- t(as.matrix(X)) 
-      
-      if(dim(X)[2] != M) {
-        stop(paste("matrix of evaluation points should be a 2 columns matrix"))
-      }
-      evals <- NULL  
-      if(length(FunctionSpace$mesh$times) == 0){
-        evals <- apply(FunctionSpace$basis_function$eval(as.matrix(X)) %*% coeff, 
-                       MARGIN=1, FUN=sum)
-      }else{
-        evals <- matrix(nrow=nrow(X),ncol=length(FunctionSpace$mesh$times))
-        for(t in 1:length(FunctionSpace$mesh$times)){
-          evals[,t] <- apply(FunctionSpace$basis_function$eval(as.matrix(X)) %*% coeff[,t], 
-                             MARGIN=1, FUN=sum) 
-        }  
-      }
-      return(evals)
-    },
-    set_coeff = function(coefficients){
-      coeff <<- coefficients
-    }
-  )
+.FunctionCtr <- R6Class("Function",
+                        public = list(
+                          FunctionSpace  = "FunctionSpace", 
+                          coeff = "matrix",
+                          initialize = function(FunctionSpace, coeff){
+                            self$FunctionSpace <- FunctionSpace
+                            self$coeff <- coeff
+                          },
+                          eval = function(X) {
+                            M = dim(self$FunctionSpace$mesh$get_nodes())[2]
+                            if(is.vector(X)) X <- t(as.matrix(X)) 
+                            
+                            if(dim(X)[2] != M) {
+                              stop(paste("matrix of evaluation points should be a 2 columns matrix"))
+                            }
+                            evals <- NULL  
+                            if(length(self$FunctionSpace$mesh$times) == 0){
+                              evals <- apply(self$FunctionSpace$basis_function$eval(as.matrix(X)) %*% self$coeff, 
+                                             MARGIN=1, FUN=sum)
+                            }else{
+                              evals <- matrix(nrow=nrow(X),ncol=length(self$FunctionSpace$mesh$times))
+                              for(t in 1:length(self$FunctionSpace$mesh$times)){
+                                evals[,t] <- apply(self$FunctionSpace$basis_function$eval(as.matrix(X)) %*% self$coeff[,t], 
+                                                   MARGIN=1, FUN=sum) 
+                              }  
+                            }
+                            return(evals)
+                          },
+                          set_coeff = function(coefficients){
+                            self$coeff <- coefficients
+                          }
+                        )
 )
 
+#' Finite Element Function
+#'
+#' @name Function
+#'
+#' @exportClass Function
+setOldClass(c("Function", "R6"))
+
 ## constructor
 
 #' Create Function object
@@ -137,7 +165,7 @@ setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
 #' }
 Function <- function(FunctionSpace) {
   coeff = matrix(ncol = 1, nrow = 0)
-  .FunctionCtr(coeff = coeff, FunctionSpace = FunctionSpace)
+  .FunctionCtr$new(coeff = coeff, FunctionSpace = FunctionSpace)
 }
 
 ## Function plot overload
diff --git a/R/mesh.R b/R/mesh.R
index 2084944..cb94265 100644
--- a/R/mesh.R
+++ b/R/mesh.R
@@ -1,37 +1,40 @@
 
 ## Mesh Class Definition
-.MeshCtr <- setRefClass(
-  Class = "Mesh", contains = "Domain",
-  fields = c(
-    cpp_handler = "ANY",       ## cpp backend
-    m = "integer",             ## local dim
-    n = "integer",             ## embedding dim
-    times = "numeric",
-    deltaT = "numeric" # for parabolic problems
-  ),
-  methods = c(
-    get_nodes = function(){
-      cpp_handler$nodes()
-    },
-    get_elements = function(){
-      cpp_handler$elements()+1
-    },
-    get_boundary = function(){
-      cpp_handler$boundary()
-    },
-    get_neighbors = function(){
-      cpp_handler$neighbors()
-    },
-    get_times = function(){
-      return(times)
-    },
-    set_deltaT = function(deltaT){
-      if(length(time_interval)==0)
-        stop("Error!")
-      deltaT <<- deltaT
-      times <<- seq(time_interval[1], time_interval[2], by=deltaT)
-    }
-  )
+#' @export 
+.MeshCtr <- R6Class("Mesh", inherit = .DomainCtr,
+                    public = list(
+                      cpp_handler = "ANY",       ## cpp backend
+                      m = 0L,             ## local dim
+                      n = 0L,             ## embedding dim
+                      times = vector(mode="double", length = 0L),
+                      deltaT = NA, # for parabolic problems
+                      initialize = function(cpp_handler, m, n){
+                        self$cpp_handler <- cpp_handler
+                        self$m <- m
+                        self$n <- n
+                      },
+                      get_nodes = function(){
+                        self$cpp_handler$nodes()
+                      },
+                      get_elements = function(){
+                        self$cpp_handler$elements()+1
+                      },
+                      get_boundary = function(){
+                        self$cpp_handler$boundary()
+                      },
+                      get_neighbors = function(){
+                        self$cpp_handler$neighbors()
+                      },
+                      get_times = function(){
+                        return(self$times)
+                      },
+                      set_deltaT = function(deltaT){
+                        if(length(self$time_interval)==0)
+                          stop("Error!")
+                        self$deltaT <- deltaT
+                        self$times <- seq(self$time_interval[1], self$time_interval[2], by=deltaT)
+                      }
+                    )
 )
 
 #' Create mesh object
@@ -67,13 +70,12 @@ setMethod("Mesh", signature = c(domain="list"),
             m <- ncol(domain$elements) - 1
             n <- ncol(domain$nodes) 
             if(m == 2 & n == 2)
-              .MeshCtr(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n), 
-                       times=vector(mode="double"))                                                  
+              .MeshCtr$new(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n))                                                  
             else
               stop("wrong input argument provided.")
 })
 
-# @importClassesFrom RTriangle triangulation
+setOldClass("triangulation")
 
 #' @rdname Mesh
 setMethod("Mesh", signature=c(domain="triangulation"),
@@ -91,9 +93,7 @@ setMethod("Mesh", signature=c(domain="triangulation"),
             
             domain <- list(elements = elements, nodes = nodes, boundary = boundary)
             if(m == 2 & n == 2)
-              .MeshCtr(cpp_handler=new(Mesh_2D, domain), m=as.integer(m),n=as.integer(n), 
-                       times=vector(mode="double"),
-                       time_interval = vector(mode="numeric", length = 0), crs = NA_crs_)                                                  
+              .MeshCtr$new(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n))                                                  
             else
               stop("wrong input argument provided.")
 })
diff --git a/R/operators.R b/R/operators.R
index 85a2964..82059bf 100644
--- a/R/operators.R
+++ b/R/operators.R
@@ -1,12 +1,22 @@
 # gradient of Function
-.FunctionGradCtr <- setRefClass(
-  Class = "FunctionGrad",
-  fields = c(
+.FunctionGradCtr <- R6Class("FunctionGrad",
+  public = list(
     f = "Function", ## Function of which the gradient is taken
-    K = "ANY" ## set by product operator
+    K = "ANY",      ## set by product operator
+    initialize = function(f, K=matrix(c(1,0,0,1), nrow=2,ncol=2,byrow=T)){
+      self$f <- f
+      self$K <- K
+    }
   )
 )
 
+#' Function gradient
+#'
+#' @name grad
+#'
+#' @exportClass FunctionGrad
+setOldClass(c("FunctionGrad", "R6"))
+
 ## take gradient of Function
 
 #' compute gradient of Function
@@ -28,31 +38,30 @@ setGeneric("grad", function(f) standardGeneric("grad"))
 
 #' @rdname grad
 setMethod("grad", signature(f = "Function"), function(f) {
-  .FunctionGradCtr(f = f)
+  .FunctionGradCtr$new(f = f)
 })
 
+# setMethod("*", signature=c(e1="matrix", e2="FunctionGrad")
+
 #' product overload for FunctionGradObejct
 #'
 #' @param e1 a numeric matrix.
-#' @param e2 a FunctionGrad created by \code{grad} function.
-#' @return A FunctionGrad.
+#' @param e2 a FunctionGrad class created by \code{grad} function.
+#' @return A FunctionGrad class.
 #' @export 
-setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"),
-          function(e1,e2){
-        .FunctionGradCtr(f = e2$f, K = e1)    
-})
-
-# ## diffusion tensor - FunctionGrad product overload
-# `*.FunctionGrad` <- function(op1, op2) {
-#   if (!is.matrix(op1) && !is.function(op1))
-#     stop("bad diffusion tensor type")
-#   .FunctionGradCtr(f = op2$f, K = op1)
-# }
+`*.FunctionGrad` <- function(e1, e2){
+  if (!is.matrix(e1)) stop("First argument must be a matrix!")
+        .FunctionGradCtr$new(f = e2$f, K = e1)    
+}
 
 ## base class for differential operators
-.DiffOpCtr <- setRefClass(
-    Class = "DiffOpObject",
-    fields = c(
+.DiffOpCtr <- R6Class("DiffOpObject",
+    public = list(
+        initialize = function(tokens, params, f){
+          self$tokens <- tokens
+          self$params <- params
+          self$f <- f
+        },
         ## single blocks composing the overall operator
         tokens = "vector",
         params = "list",
@@ -61,6 +70,15 @@ setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"),
     )
 )
 
+#' R6 class representing a differential operator
+#' 
+#' @name Differential Operator
+#'
+#' @exportClass DiffOpObject
+setOldClass(c("DiffOpObject", "R6"))
+
+##setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
+
 ## sum of differential operators
 
 #' plus operator overload for DiffOpObject
@@ -68,23 +86,24 @@ setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"),
 #' @param e1 a DiffOpObject.
 #' @param e2 a DiffOpObject.
 #' @return A S4 object representing the sum of two differential operators.
+#' @name DifferentialOps
 #' @export 
-setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
-          function(e1, e2) {
-            if (tracemem(e1$f) != tracemem(e2$f)) {
+`+.DiffOpObject` <- function(e1, e2){
+            if (!identical(e1$f, e2$f)) {
               stop("operator arguments must be the same")
             }
             ## check for duplicated operator tokens
             if (any(duplicated(c(e1$tokens, e2$tokens)))) {
               stop("duplicated operator")
             }
-            .DiffOpCtr(
+            .DiffOpCtr$new(
               tokens = c(e1$tokens, e2$tokens),
               params = c(e1$params, e2$params),
               f = e1$f
             )
-          }
-)
+}
+
+# setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
 
 #' difference of differential operators
 #'
@@ -92,10 +111,13 @@ setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
 #' @param e2 a DiffOpObject.
 #' @return A S4 object representing the difference of two differential operators.
 #' @rdname minus_DiffOb_op
+#' @name DifferentialOps
 #' @export 
-setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
-          function(e1, e2) {
-            if (tracemem(e1$f) != tracemem(e2$f)) {
+`-.DiffOpObject` <- function(e1, e2){
+            if(missing(e2)){
+              e1$params[[1]] <- -e1$params[[1]]
+              return(e1)
+            }else if (!identical(e1$f, e2$f)) {
               stop("operator arguments must be the same")
             }
             ## check for duplicated operator tokens
@@ -105,55 +127,62 @@ setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
             for(i in 1:length(e2$params)){
               e2$params[[i]] = -e2$params[[i]] 
             }
-            .DiffOpCtr(
+            .DiffOpCtr$new(
               tokens = c(e1$tokens, e2$tokens),
               params = c(e1$params, e2$params),
               f = e1$f
             )
-})
+  }
+
+#setMethod("-", signature =c(e1 = "DiffOpObject", e2 = "missing"),
 
- 
-# minus (unary) operator for DiffOpObject
+# minus (unary) operator for differential operators
 # 
 # @param e1 a DiffOpObject.
 # @return A S4 object differential operator whose paramter has changed sign.
 # @export
 
-#' @rdname minus_DiffOb_op
-setMethod("-", signature =c(e1 = "DiffOpObject", e2 = "missing"),
-          function(e1, e2){
-            e1$params[[1]] <- -e1$params[[1]]
-            e1
-          }
-)
+# @rdname minus_DiffOb_op
+# `-.DiffOpObject` <- function(e1){
+#             e1$params[[1]] <- -e1$params[[1]]
+#             e1
+# }
+
+# setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"),
 
 #' 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
 #' @export 
-setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"),
-          function(e1,e2){
-            #if (!is.numeric(e1)) stop("bad product")
+`*.DiffOpObject` <- function(e1,e2){  
+            if (!is.numeric(e1)) stop("bad product")
             e2$params[[1]] <- e1*e2$params[[1]]
             e2
-            
-})
-
+          }
+  
 ## diffusion term
-.DiffusionCtr <- setRefClass(
-    Class = "DiffusionOperator",
-    contains = "DiffOpObject"
+.DiffusionCtr <- R6Class("DiffusionOperator",
+    inherit = .DiffOpCtr
 )
 
+#' Diffusion Operator
+#'
+#' @name DifferentialOperator
+#'
+#' @exportClass DiffusionOperator
+setOldClass(c("DiffusionOperator", "DiffOpObject"))
+
 ## laplace() returns a special operator for the case of
 ## isotropic  and stationary diffusion
 
-#' laplace operator for Function
+#' laplace operator for Function class
 #'
 #' @param f a Function.
 #' @return A S4 object representing the laplace operator applied to the function passed as parameter.
+#' @rdname DiffusionOperator
 #' @export 
 #' @examples
 #' \dontrun{
@@ -168,7 +197,7 @@ laplace <- function(f) {
     if (!is(f, "Function")) {
         stop("wrong argument type")
     }
-    .DiffusionCtr(
+    .DiffusionCtr$new(
         tokens = "diffusion",
         params = list(diffusion = 1),
         f = f
@@ -181,6 +210,7 @@ laplace <- function(f) {
 #'
 #' @param f a Function.
 #' @return A S4 object representing the diffusion term of a second order linear differential operator.
+#' @rdname DiffusionOperator
 #' @export
 #' @examples
 #' \dontrun{
@@ -195,7 +225,7 @@ laplace <- function(f) {
 div <- function(f) {
     if (is(f, "FunctionGrad")) {
         if (!is.null(f$K)) {
-            return(.DiffusionCtr(
+            return(.DiffusionCtr$new(
                        tokens = "diffusion",
                        params = list(diffusion = f$K),
                        f = f$f)
@@ -206,17 +236,23 @@ div <- function(f) {
 }
 
 ## transport term
-.TransportCtr <- setRefClass(
-    Class = "TransportOperator",
-    contains = "DiffOpObject"
+.TransportCtr <- R6Class("TransportOperator",
+    inherit = .DiffOpCtr
 )
 
+#' Transport Operator
+#'
+#' @name TransportOperator
+#'
+#' @exportClass TransportOperator
+setOldClass(c("TransportOperator", "DiffOpObject"))
+
 #' 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.
-#' @rdname dot_product
+#' @rdname TransportOperator
 #' @export
 #' @examples
 #' \dontrun{
@@ -230,10 +266,10 @@ div <- function(f) {
 #' } 
 setGeneric("dot", function(op1, op2) standardGeneric("dot"))
 
-#' @rdname dot_product
+#' @rdname TransportOperator
 setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"),
           function(op1, op2) {
-              .TransportCtr(
+              .TransportCtr$new(
                   tokens = "transport",
                   params = list(transport = as.matrix(op1)),
                   f = op2$f
@@ -241,16 +277,25 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"),
 })
 
 ## reaction term
-.ReactionCtr <- setRefClass(
-    Class = "ReactionOperator",
-    contains = "DiffOpObject"
+.ReactionCtr <- R6Class("ReactionOperator",
+    inherit = .DiffOpCtr
 )
 
+#' Reaction Operator
+#'
+#' @name ReactionOperator
+#'
+#' @exportClass ReactionOperator
+setOldClass(c("ReactionOperator", "DiffOpObject"))
+
+# setMethod("*", signature = c(e1="numeric", e2="Function"),
+          
 #' 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.
+#' @rdname ReactionOperator
 #' @export 
 #' @examples
 #' \dontrun{
@@ -261,20 +306,24 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"),
 #' f <- Function(Vh)
 #' reaction <- 2*f
 #' }
-setMethod("*", signature = c(e1="numeric", e2="Function"),
-          function(e1,e2){
-            .ReactionCtr(
+`*.Function` <- function(e1,e2){
+            if(!is(e1,"numeric")) stop("First argument must be a scalar!")
+            .ReactionCtr$new(
               tokens = "reaction",
               params = list(reaction = e1),
               f = e2
             )
-})
+}
 
-.TimeDerivativeCtr <- setRefClass(
-  Class = "TimeDerivative",
-  contains ="DiffOpObject"
+.TimeDerivativeCtr <- R6Class("TimeDerivative",
+  inherit = .DiffOpCtr
 )
 
+#' @name TimeDerivative
+#'
+#' @exportClass TimeDerivative
+setOldClass(c("TimeDerivative", "DiffOpObject"))
+
 # overloading stats::dt function
 
 #' time derivate of FunctionObejct
@@ -283,6 +332,7 @@ setMethod("*", signature = c(e1="numeric", e2="Function"),
 #' @param df missing.
 #' @param ncp missing.
 #' @return A S4 object representing the time derivative of a Function.
+#' @rdname TimeDerivative
 #' @export
 #' @examples
 #' \dontrun{
@@ -295,7 +345,7 @@ setMethod("*", signature = c(e1="numeric", e2="Function"),
 #' }
 setMethod("dt", signature = c(x="Function", df="missing", ncp="missing"),
            function(x,df,ncp){
-             .TimeDerivativeCtr(
+             .TimeDerivativeCtr$new(
                tokens="time",
                params = list(time=1L),
               f=x
diff --git a/R/pde.R b/R/pde.R
index 877f7bd..178096c 100644
--- a/R/pde.R
+++ b/R/pde.R
@@ -2,123 +2,141 @@
 pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3)
 
 # Pde Class wraps C++ R_PDE class
-.PdeCtr <- setRefClass(
-  Class = "PdeObject",
-  fields = c(
+.PdeCtr <- R6Class("Pde",
+  public = list(
     is_dirichletBC_set = "logical",      
     is_initialCondition_set = "logical",
     is_parabolic ="logical",          
-    cpp_handler = "ANY"),
-  methods = c(
+    cpp_handler = "ANY",
+    DifferentialOp = "DiffOpObject",
+    initialize = function(cpp_handler, DifferentialOp, is_parabolic, is_dirichletBC_set, is_initialCondition_set){
+      self$cpp_handler <- cpp_handler
+      self$DifferentialOp <- DifferentialOp 
+      self$is_parabolic <- is_parabolic
+      self$is_dirichletBC_set <- is_dirichletBC_set
+      self$is_initialCondition_set <- is_initialCondition_set
+    },
     solve = function(){
-      if(!is_dirichletBC_set){
-        if(!is_parabolic){
-          dirichletBC_ = as.matrix(rep(0,times=nrow(cpp_handler$get_dofs_coordinates())))
+      if(!self$is_dirichletBC_set){
+        if(!self$is_parabolic){
+          dirichletBC_ = as.matrix(rep(0,times=nrow(self$cpp_handler$get_dofs_coordinates())))
         }else 
-          dirichletBC_ = matrix(0, nrow=nrow(cpp_handler$get_dofs_coordinates()), ncol=length(times))
-      cpp_handler$set_dirichlet_bc(dirichletBC_)
-      is_dirichletBC_set <<- TRUE
+          dirichletBC_ = matrix(0, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), 
+                                   ncol=length(self$DifferentialOp$f$FunctionSpace$mesh$times))
+      self$cpp_handler$set_dirichlet_bc(dirichletBC_)
+      self$is_dirichletBC_set <- TRUE
       }
-      if(is_parabolic & (!is_initialCondition_set)){
+      if(self$is_parabolic & (!self$is_initialCondition_set)){
         stop("initialCondition must be provided.")
       }
-      cpp_handler$solve()
-    },
-    get_solution = function(){
-      cpp_handler$solution()
+      self$cpp_handler$solve()
     },
     get_dofs_coordinates = function(){
-      cpp_handler$get_dofs_coordinates()
+      self$cpp_handler$get_dofs_coordinates()
     },
     set_boundary_condition = function(fun, type="dirichlet", on=NULL){
       if(!any(type == c("dirichlet", "Dirichlet", "d"))) 
         stop("Only Dirichlet boundary condtions allowed.")
       dirichletBC_ <- NULL
-      if(typeof(fun) == "closure"){
-        if(!is_parabolic){
-          dirichletBC_ <- as.matrix(fun(cpp_handler$get_dofs_coordinates()))
+      if(any(typeof(fun) == c("function", "closure"))){
+        if(!self$is_parabolic){
+          dirichletBC_ <- as.matrix(fun(self$cpp_handler$get_dofs_coordinates()))
         }else{
-          dirichletBC_ <- fun(cpp_handler$get_dofs_coordinates(),times)
+          dirichletBC_ <- fun(self$cpp_handler$get_dofs_coordinates(), 
+                              self$DifferentialOp$f$FunctionSpace$mesh$times)
         }
       }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){
-        if(nrow(as.matrix(fun)) == nrow(cpp_handler$get_dofs_coordinates())){
+        if(nrow(as.matrix(fun)) == nrow(self$cpp_handler$get_dofs_coordinates())){
           dirichletBC_ <- fun  
         }else if (nrow(as.matrix(fun)) == 1L){
-          if(!is_parabolic)
-            dirichletBC_ <- matrix(fun, nrow=nrow(cpp_handler$get_dofs_coordinates()), ncol=1)
+          if(!self$is_parabolic)
+            dirichletBC_ <- matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), ncol=1)
           else 
-            dirichletBC_ <- matrix(fun, nrow=nrow(cpp_handler$get_dofs_coordinates()), ncol=length(times))
+            dirichletBC_ <- matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), 
+                                        ncol=length(self$DifferentialOp$f$FunctionSpace$mesh$times))
           }   
         }
-      is_dirichletBC_set <<- TRUE
-      cpp_handler$set_dirichlet_bc(dirichletBC_)
+      self$is_dirichletBC_set <- TRUE
+      self$cpp_handler$set_dirichlet_bc(dirichletBC_)
     },
     set_initial_condition = function(fun){
-      if(!is_parabolic)
+      if(!self$is_parabolic)
         stop("Cannot set initial condition for elliptic problem.")
-      is_initialCondition_set <<- TRUE
+      self$is_initialCondition_set <- TRUE
       if(typeof(fun) == "closure"){ 
-        cpp_handler$set_initial_condition(fun(cpp_handler$get_dofs_coordinates()))
+        self$cpp_handler$set_initial_condition(fun(self$cpp_handler$get_dofs_coordinates()))
       }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){
-        if(nrow(as.matrix(fun)) == nrow(cpp_handler$get_dofs_coordinates())){
-          cpp_handler$set_initial_condition(fun)
+        if(nrow(as.matrix(fun)) == nrow(self$cpp_handler$get_dofs_coordinates())){
+          self$cpp_handler$set_initial_condition(fun)
         }else if(nrow(as.matrix(fun)) == 1L){
-          cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(cpp_handler$get_dofs_coordinates()),
+          self$cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()),
                                                  ncol=1))   
         }
       }
     },
     get_mass = function(){
-      cpp_handler$mass()
+      self$cpp_handler$mass()
     },
     get_stiff = function(){
-      cpp_handler$stiff()
+      self$cpp_handler$stiff()
     }
   )
 )
 
+#' @name pde
+#'
+#' @exportClass Pde
+setOldClass(c("Pde", "R6"))
+
 ## infers the type of a pde
-extract_pde_type <- function(L) {
-  if ("time" %in% names(L$params)) {
+extract_pde_type <- function(DifferentialOp) {
+  if ("time" %in% names(DifferentialOp$params)) {
     return(pde_type_list$parabolic)
   }
-  if ("diffusion" %in% names(L$params) && !is.matrix(L$params$diffusion)) {
+  if ("diffusion" %in% names(DifferentialOp$params) && !is.matrix(DifferentialOp$params$diffusion)) {
     return(pde_type_list$laplacian)
   }
   return(pde_type_list$elliptic)
 }
 
 ## parse pde parameters
-parse_pde_parameters <- function(L){
+parse_pde_parameters <- function(DifferentialOp){
+  pde_type <- extract_pde_type(DifferentialOp)
+  
   pde_parameters <- NULL
   pde_parameters$diffusion <- 1.0
   pde_parameters$transport <- matrix(0, nrow = 2, ncol = 1)
   pde_parameters$reaction  <- 0.0
   
-  for (i in length(L$params)) {
-    pde_parameters[[paste(L$tokens[i])]] <- L$params[[paste(L$tokens[i])]]
+  for (i in 1:length(DifferentialOp$params)) {
+    pde_parameters[[DifferentialOp$tokens[i]]] <- DifferentialOp$params[[DifferentialOp$tokens[i]]]
   }
+  
+  if(pde_type == pde_type_list$parabolic)
+    pde_parameters[["time_mesh"]] <- as.vector(DifferentialOp$f$FunctionSpace$mesh$times)
+  
+  if(pde_type == pde_type_list$parabolic & is(pde_parameters[["diffusion"]], "numeric"))
+    pde_parameters[["diffusion"]] <- pde_parameters[["diffusion"]]*matrix(c(1,0,0,1), nrow=2, ncol=2, byrow=T)
+  
   return(pde_parameters)
 }
 
 ## build cpp object
-make_pde <- function(L) {
-  pde_type = extract_pde_type(L)
+make_pde <- function(DifferentialOp) {
+  pde_type <- extract_pde_type(DifferentialOp)
   
   ## pde parameters
-  pde_parameters <- parse_pde_parameters(L)
-  
-  if(pde_type == pde_type_list$parabolic) pde_parameters["time_mesh"] <- L$f$FunctionSpace$mesh$times
+  pde_parameters <- parse_pde_parameters(DifferentialOp)
   
   ## define Rcpp module
-  D <- L$f$FunctionSpace$mesh$cpp_handler ## domain
-  fe_order <- L$f$FunctionSpace$fe_order  ## finite element order
+  D <- DifferentialOp$f$FunctionSpace$mesh$cpp_handler ## domain
+  fe_order <- DifferentialOp$f$FunctionSpace$fe_order  ## finite element order
   cpp_handler <- NULL
   if (fe_order == 1) { ## linear finite elements
-    cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters, L$f)
+    cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters, DifferentialOp$f)
   }
   if (fe_order == 2) { ## quadratic finite elements
-    cpp_handler <- new(cpp_pde_2d_fe2, D, pde_type - 1, pde_parameters, L$f)
+    cpp_handler <- new(cpp_pde_2d_fe2, D, pde_type - 1, pde_parameters, DifferentialOp$f)
   }
   
   return(cpp_handler)
@@ -126,24 +144,24 @@ make_pde <- function(L) {
 
 #' A PDEs object
 #'
-#' @param L a differential operator.
+#' @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.
-#' @return A S4 object representing a PDE.
+#' @return A R6 class representing a PDE.
 #' @rdname pde
 #' @export 
-setGeneric("Pde", function(L,f) standardGeneric("Pde"))
+setGeneric("Pde", function(DifferentialOp,f) standardGeneric("Pde"))
 
 #' @rdname pde
-setMethod("Pde", signature=c(L="DiffOpObject", f="ANY"),
-          function(L,f){
+setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"),
+          function(DifferentialOp,f){
             
-            pde_type = extract_pde_type(L)
-            cpp_handler <- make_pde(L)
+            pde_type = extract_pde_type(DifferentialOp)
+            cpp_handler <- make_pde(DifferentialOp)
             
             quad_nodes <- cpp_handler$get_quadrature_nodes()
             ## evaluate forcing term on quadrature nodes
             if(pde_type == pde_type_list$parabolic) {
-              times <- L$f$FunctionSpace$mesh$times
+              times <- DifferentialOp$f$FunctionSpace$mesh$times
               cpp_handler$set_forcing(as.matrix(f(quad_nodes, times)))
             }else{
               cpp_handler$set_forcing(as.matrix(f(quad_nodes)))
@@ -155,16 +173,17 @@ setMethod("Pde", signature=c(L="DiffOpObject", f="ANY"),
             is_parabolic = pde_type == pde_type_list$parabolic
             
             ## return
-            .PdeCtr(is_dirichletBC_set = FALSE,
-                    is_initialCondition_set = FALSE,
-                    is_parabolic = is_parabolic,
-                    cpp_handler = cpp_handler)
+            .PdeCtr$new(is_dirichletBC_set = FALSE,
+                        is_initialCondition_set = FALSE,
+                        is_parabolic = is_parabolic,
+                        cpp_handler = cpp_handler,
+                        DifferentialOp = DifferentialOp)
 })
 
 #' @rdname pde
-setMethod("Pde", signature=c(L="DiffOpObject", f="numeric"),
-          function(L,f){
-            pde_type = extract_pde_type(L)
+setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="numeric"),
+          function(DifferentialOp,f){
+            pde_type = extract_pde_type(DifferentialOp)
             
             fun <- NULL
             if(pde_type!=pde_type_list$parabolic){
@@ -177,5 +196,5 @@ setMethod("Pde", signature=c(L="DiffOpObject", f="numeric"),
               }  
             }
             
-            return(Pde(L,fun))
+            return(Pde(DifferentialOp,fun))
 })
diff --git a/man/div.Rd b/man/div.Rd
deleted file mode 100644
index d679390..0000000
--- a/man/div.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{div}
-\alias{div}
-\title{divergence operator FunctionGrad}
-\usage{
-div(f)
-}
-\arguments{
-\item{f}{a Function.}
-}
-\value{
-A S4 object representing the diffusion term of a second order linear differential operator.
-}
-\description{
-divergence operator FunctionGrad
-}
-\examples{
-\dontrun{
-library(femR)
-data("unit_square")
-mesh <- Mesh(unit_square)
-Vh <- FunctionSpace(mesh)
-f <- Function(Vh)
-K <- matrix(c(1,2,1,0),nrow=2,ncol=2)
-diffusion <- div(K*grad(f))
-} 
-}
diff --git a/man/dot_product.Rd b/man/dot_product.Rd
deleted file mode 100644
index 30d7bb7..0000000
--- a/man/dot_product.Rd
+++ /dev/null
@@ -1,33 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{dot}
-\alias{dot}
-\alias{dot,vector,FunctionGrad-method}
-\title{dot product between vector and FunctionGrad}
-\usage{
-dot(op1, op2)
-
-\S4method{dot}{vector,FunctionGrad}(op1, op2)
-}
-\arguments{
-\item{op1}{a numeric vector.}
-
-\item{op2}{a FunctionGrad.}
-}
-\value{
-A S4 object representing the advection term of a second order linear differential operator.
-}
-\description{
-dot product between vector and FunctionGrad
-}
-\examples{
-\dontrun{
-library(femR)
-data("unit_square")
-mesh <- Mesh(unit_square)
-Vh <- FunctionSpace(mesh)
-f <- Function(Vh)
-b <- c(1,1)
-advection <- dot(b,grad(f))
-} 
-}
diff --git a/man/dt-Function-missing-missing-method.Rd b/man/dt-Function-missing-missing-method.Rd
deleted file mode 100644
index f9be2af..0000000
--- a/man/dt-Function-missing-missing-method.Rd
+++ /dev/null
@@ -1,31 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{dt,Function,missing,missing-method}
-\alias{dt,Function,missing,missing-method}
-\title{time derivate of FunctionObejct}
-\usage{
-\S4method{dt}{Function,missing,missing}(x, df, ncp)
-}
-\arguments{
-\item{x}{a Function.}
-
-\item{df}{missing.}
-
-\item{ncp}{missing.}
-}
-\value{
-A S4 object representing the time derivative of a Function.
-}
-\description{
-time derivate of FunctionObejct
-}
-\examples{
-\dontrun{
-library(femR)
-data("unit_square")
-mesh <- Mesh(unit_square)
-Vh <- FunctionSpace(mesh)
-f <- Function(Vh)
-dt(f)
-}
-}
diff --git a/man/laplace.Rd b/man/laplace.Rd
deleted file mode 100644
index 5c6e12f..0000000
--- a/man/laplace.Rd
+++ /dev/null
@@ -1,27 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{laplace}
-\alias{laplace}
-\title{laplace operator for Function}
-\usage{
-laplace(f)
-}
-\arguments{
-\item{f}{a Function.}
-}
-\value{
-A S4 object representing the laplace operator applied to the function passed as parameter.
-}
-\description{
-laplace operator for Function
-}
-\examples{
-\dontrun{
-library(femR)
-data("unit_square")
-mesh <- Mesh(unit_square)
-Vh <- FunctionSpace(mesh)
-f <- Function(Vh)
-laplace_f <- laplace(f)
-}
-}
diff --git a/man/plus-DiffOpObject-DiffOpObject-method.Rd b/man/plus-DiffOpObject-DiffOpObject-method.Rd
deleted file mode 100644
index 659e296..0000000
--- a/man/plus-DiffOpObject-DiffOpObject-method.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{+,DiffOpObject,DiffOpObject-method}
-\alias{+,DiffOpObject,DiffOpObject-method}
-\title{plus operator overload for DiffOpObject}
-\usage{
-\S4method{+}{DiffOpObject,DiffOpObject}(e1, e2)
-}
-\arguments{
-\item{e1}{a DiffOpObject.}
-
-\item{e2}{a DiffOpObject.}
-}
-\value{
-A S4 object representing the sum of two differential operators.
-}
-\description{
-plus operator overload for DiffOpObject
-}
diff --git a/man/times-matrix-FunctionGrad-method.Rd b/man/times-matrix-FunctionGrad-method.Rd
deleted file mode 100644
index 1033a27..0000000
--- a/man/times-matrix-FunctionGrad-method.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{*,matrix,FunctionGrad-method}
-\alias{*,matrix,FunctionGrad-method}
-\title{product overload for FunctionGradObejct}
-\usage{
-\S4method{*}{matrix,FunctionGrad}(e1, e2)
-}
-\arguments{
-\item{e1}{a numeric matrix.}
-
-\item{e2}{a FunctionGrad created by \code{grad} function.}
-}
-\value{
-A FunctionGrad.
-}
-\description{
-product overload for FunctionGradObejct
-}
diff --git a/man/times-numeric-DiffOpObject-method.Rd b/man/times-numeric-DiffOpObject-method.Rd
deleted file mode 100644
index 9101d62..0000000
--- a/man/times-numeric-DiffOpObject-method.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{*,numeric,DiffOpObject-method}
-\alias{*,numeric,DiffOpObject-method}
-\title{product by scalar for differential operators}
-\usage{
-\S4method{*}{numeric,DiffOpObject}(e1, e2)
-}
-\arguments{
-\item{e1}{a numeric.}
-
-\item{e2}{a DiffOpObject.}
-}
-\value{
-A S4 object representing the product by scalar for a differential operator.
-}
-\description{
-product by scalar for differential operators
-}
diff --git a/man/times-numeric-Function-method.Rd b/man/times-numeric-Function-method.Rd
deleted file mode 100644
index c56f218..0000000
--- a/man/times-numeric-Function-method.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/operators.R
-\name{*,numeric,Function-method}
-\alias{*,numeric,Function-method}
-\title{product by scalar for FunctionObejct}
-\usage{
-\S4method{*}{numeric,Function}(e1, e2)
-}
-\arguments{
-\item{e1}{a numeric.}
-
-\item{e2}{a FunctioObject created by \code{Function}.}
-}
-\value{
-A S4 object representing the reaction term of a second order linear differential operator.
-}
-\description{
-product by scalar for FunctionObejct
-}
-\examples{
-\dontrun{
-library(femR)
-data("unit_square")
-mesh <- Mesh(unit_square)
-Vh <- FunctionSpace(mesh)
-f <- Function(Vh)
-reaction <- 2*f
-}
-}
diff --git a/tests/adv.diff.K.R b/tests/adv.diff.K.R
index ee72fe1..8fabddd 100644
--- a/tests/adv.diff.K.R
+++ b/tests/adv.diff.K.R
@@ -34,7 +34,7 @@ f <- Function(Vh)
 #L <- -1*laplace(f) + dot(c(alpha_,0), grad(f)) # OK
 I <- matrix(c(1,0,0,1), byrow=T, nrow=2,ncol=2)
 L <- -div(I*grad(f)) + dot(c(alpha_,0), grad(f))
-
+# L <- -div(grad(f)) + dot(c(alpha_,0), grad(f))
 ## forcing term
 u <- function(points){
   return(gamma_ * sin(pi * points[,2])) 
@@ -53,18 +53,18 @@ pde$solve()
 
 ## compute L2 norm of the error
 u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates()))
-error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2))
+error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - f$coeff)^2))
 cat("L2 error = ", error.L2, "\n")
 
 ## perform evaluation at single point
 point = c(0.2, 0.5)
-f$eval_at(point)
+f$eval(point)
 
 ## evaluate over a 10x10 grid
 x <- seq(0, 1, length.out = 50)
 y <- x
 points <- expand.grid(x, y)
-max(abs( f$eval_at(points) - as.matrix(exact_solution(points))))
+max(abs( f$eval(points) - as.matrix(exact_solution(points))))
 
 # plot solution 
 options(warn=-1)
diff --git a/tests/adv.diff.R b/tests/adv.diff.R
index 717243b..e23d5e6 100644
--- a/tests/adv.diff.R
+++ b/tests/adv.diff.R
@@ -47,19 +47,19 @@ pde$solve()
 
 ## compute L2 norm of the error
 u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates()))
-error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2))
+error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - f$coeff)^2))
 cat("L2 error = ", error.L2, "\n")
 
 ## perform evaluation at single point
 point = c(0.2, 0.5)
-f$eval_at(point)
+f$eval(point)
 
 ## evaluate over a 10x10 grid
 x <- seq(0, 1, length.out = 50)
 y <- x
 points <- expand.grid(x, y)
-f$eval_at(points)
-max(abs( f$eval_at(points) - as.matrix(exact_solution(points))))
+f$eval(points)
+max(abs( f$eval(points) - as.matrix(exact_solution(points))))
 # plot solution 
 options(warn=-1)
 
diff --git a/tests/parabolic.R b/tests/parabolic.R
index d0fc942..6ce7336 100644
--- a/tests/parabolic.R
+++ b/tests/parabolic.R
@@ -33,7 +33,7 @@ u <- Function(Vh)
 ## define differential operator in its strong formulation
 #L <- dt(f) + (-1)*laplace(f) 
 
-L <- dt(u) - laplace(u) + dot(c(0,10),grad(u))
+L <- dt(u) - laplace(u) + dot(c(0,1),grad(u))
 
 ## forcing term
 f <- function(points,times){
@@ -60,6 +60,7 @@ pde$set_initial_condition(u0)
 
 pde$set_boundary_condition(g(pde$get_dofs_coordinates(),times))
 pde$set_initial_condition(u0(pde$get_dofs_coordinates()))
+
 pde$solve()
 
 ## compute L2 norm of the error
@@ -67,20 +68,20 @@ u_ex <- exact_solution(pde$get_dofs_coordinates(), times)
  
 error.L2 <- matrix(0, nrow=length(times), ncol=1)
 for( t in 1:length(times)){
-  error.L2[t] <- sqrt(sum(pde$get_mass() %*% (u_ex[,t] - pde$get_solution()[,t])^2))
+  error.L2[t] <- sqrt(sum(pde$get_mass() %*% (u_ex[,t] - u$coeff[,t])^2))
   cat(paste0("L2 error at time ",times[t]," ", error.L2[t], "\n"))
 }
 # ------------------------------------------------------------------------------
 ## perform evaluation at single point
 point = c(0.2, 0.5)
-u$eval_at(point)
+u$eval(point)
 # 
 # ## evaluate over a 10x10 grid 
 x <- seq(0, 1, length.out = 50)
 y <- x
 points <- expand.grid(x, y)
-dim( u$eval_at(points) )
-max(abs( u$eval_at(points) - as.matrix(exact_solution(points))))
+dim( u$eval(points) )
+max(abs( u$eval(points) - as.matrix(exact_solution(points))))
 # 
 # ## plot solution 
 options(warn=-1)