diff --git a/.github/workflows/document.yaml b/.github/workflows/document.yaml
index 6c118dd..fec6ca6 100644
--- a/.github/workflows/document.yaml
+++ b/.github/workflows/document.yaml
@@ -15,6 +15,7 @@ jobs:
       - name: Checkout repo
         uses: actions/checkout@v3
         with:
+          submodules: 'recursive'
           fetch-depth: 0
 
       - name: Setup R
diff --git a/NAMESPACE b/NAMESPACE
index 330e930..0a56ae8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,8 +1,21 @@
 # Generated by roxygen2: do not edit by hand
 
-exportMethods(Mesh)
+export("%X%")
+export(Function)
+export(FunctionSpace)
+export(Mesh)
+export(Pde)
+export(div)
+export(dot)
+export(grad)
+export(laplace)
+export(read_mesh)
+exportMethods("*")
+exportMethods("+")
+exportMethods("-")
+exportMethods(contour)
+exportMethods(dt)
 exportMethods(plot)
-exportPattern("^[[:alpha:]]+")
 import(Rcpp)
 import(RcppEigen)
 import(methods)
diff --git a/R/function_space.R b/R/function_space.R
index faa8f0a..095249d 100644
--- a/R/function_space.R
+++ b/R/function_space.R
@@ -6,20 +6,23 @@
   )
 )
 
-# constructor
-setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpace"))
-
 #' 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.
 #' @export 
+#' @rdname FunctionSpace
+setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpace"))
+
+#' @rdname FunctionSpace
 #' @examples
+#' \dontrun{
 #' library(femR)
 #' data("unit_square")
 #' mesh <- Mesh(unit_square)
-#' Vh <- FunctionSpace(mesh)
+#' Vh <- FunctionSpace(mesh = mesh, fe_order = 1)
+#' }
 setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="numeric"),
           function(mesh,fe_order){
             if(fe_order == 1){
@@ -30,6 +33,14 @@ setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="numeric"),
                                                   
 })
 
+#' @rdname FunctionSpace
+#' @examples
+#' \dontrun{
+#' library(femR)
+#' data("unit_square")
+#' mesh <- Mesh(unit_square)
+#' Vh <- FunctionSpace(mesh)
+#' }
 setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
           function(mesh){
               return(.FunctionSpaceCtr(mesh=mesh, fe_order=1L))
@@ -66,11 +77,13 @@ setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
 #' @return An S4 object representing a Function belonging to the FunctionSpace passed as parameter.
 #' @export 
 #' @examples
+#' \dontrun{
 #' library(femR)
 #' data("unit_square")
 #' mesh <- Mesh(unit_square)
 #' Vh <- FunctionSpace(mesh)
 #' f <- Function(Vh)
+#' }
 Function <- function(FunctionSpace) {
   coeff = matrix(ncol = 1, nrow = 0)
   .FunctionCtr(coeff = coeff, FunctionSpace = FunctionSpace)
@@ -81,8 +94,8 @@ Function <- function(FunctionSpace) {
 #' Plot a Function object
 #'
 #' @param x A \code{FunctionObject} generated by \code{Function}
-#' @param ... Arguments representing graphical options to be passed to \link[plotly]{plot_ly}.
-#' @return A \link[plotly] object
+#' @param ... Arguments representing graphical options to be passed to plot_ly function.
+#' @return A plotly object
 #' @export
 setMethod("plot", signature=c(x="FunctionObject"), function(x, ...){
   times <- x$FunctionSpace$mesh$times
@@ -146,8 +159,8 @@ setMethod("plot", signature=c(x="FunctionObject"), function(x, ...){
 #' Create a contour plot of a FunctionObject
 #'
 #' @param x A \code{FunctionObject} generated by \code{Function}
-#' @param ... Arguments representing graphical options to be passed to \link[plotly]{plot_ly}.
-#' @return A \link[plotly] object
+#' @param ... Arguments representing graphical options to be passed to plot_ly function.
+#' @return A plotly object
 #' @export
 setMethod("contour", signature=c(x="FunctionObject"), function(x, ...){
   times <- x$FunctionSpace$mesh$times
diff --git a/R/mesh.R b/R/mesh.R
index b6b38af..3fa0d42 100644
--- a/R/mesh.R
+++ b/R/mesh.R
@@ -27,8 +27,6 @@
   )
 )
 
-setGeneric("Mesh", function(domain) standardGeneric("Mesh"))
-
 #' Create mesh object
 #'
 #' @param domain A named list containing must contain:
@@ -39,11 +37,11 @@ setGeneric("Mesh", function(domain) standardGeneric("Mesh"))
 #'           an entry '0' indicates that the corresponding node is not a boundary node.}
 #' }
 #' @return An S4 object representing a Mesh.
+#' @rdname MeshObject
 #' @export 
-#' @examples
-#' library(femR)
-#' data("unit_square")
-#' mesh <- Mesh(unit_square)
+setGeneric("Mesh", function(domain) standardGeneric("Mesh"))
+
+#' @rdname MeshObject
 setMethod("Mesh", signature = c(domain="list"),
           function(domain){
             domain$elements <- domain$elements - 1
@@ -57,7 +55,16 @@ setMethod("Mesh", signature = c(domain="list"),
               stop("wrong input argument provided.")
 })
 
+#' Spatio-temporal domain
+#'
+#' @param op1 A mesh object created by \code{Mesh}.
+#' @param op2 A numeric vector.
+#' @return An S4 object representing a spatio-temporal domain.
+#' @rdname MeshObject_times_vector
+#' @export 
 setGeneric("%X%", function(op1, op2) standardGeneric("%X%"))
+
+#' @rdname MeshObject_times_vector
 setMethod("%X%", signature=c(op1="MeshObject", op2="numeric"),
           function(op1, op2){
             if(op2[1] > op2[length(op2)])
@@ -119,13 +126,15 @@ plot_mesh_aux <- function(x, ...){
 #'
 #' @param x A \code{MeshObject} object defining the triangular mesh, as generated by \code{Mesh}
 #' @param ... Arguments representing graphical options to be passed to \link[plotly]{plot_ly}.
-#' @return A \link[plotly] object
+#' @return A plotly object
 #' @export
 #' @examples
+#' \dontrun{
 #' library(femR)
 #' data("unit_square")
 #' mesh <- Mesh(unit_square)
 #' plot(mesh)
+#' }
 setMethod("plot", "MeshObject", function(x, ...){
   plot_mesh_aux(x, ...)  
 })
diff --git a/R/modules.R b/R/modules.R
index 20ed7d5..26f67e9 100644
--- a/R/modules.R
+++ b/R/modules.R
@@ -2,9 +2,11 @@
 #' @import methods Rcpp
 #' @import RcppEigen
 #' @importFrom Rcpp evalCpp
-#' @exportPattern "^[[:alpha:]]+"
 NULL
 
+# #' @exportPattern "^[[:alpha:]]+"
+
+
 ## load required modules
 Rcpp::loadModule("PDE_2D_ORDER_1", TRUE)
 Rcpp::loadModule("PDE_2D_ORDER_2", TRUE)
diff --git a/R/operators.R b/R/operators.R
index 13a8c68..ce4782a 100644
--- a/R/operators.R
+++ b/R/operators.R
@@ -8,11 +8,35 @@
 )
 
 ## take gradient of Function
+
+#' Take gradient of FunctionObject
+#'
+#' @param f a FunctionObject created by \code{Function}:
+#' @return An S4 object representing the gradient of the FunctionObject passed as parameter.
+#' @export 
+#' @rdname grad
+#' @examples
+#' \dontrun{
+#' library(femR)
+#' data("unit_square")
+#' mesh <- Mesh(unit_square)
+#' Vh <- FunctionSpace(mesh)
+#' f <- Function(Vh)
+#' grad_f <- grad(f)
+#' }
 setGeneric("grad", function(f) standardGeneric("grad"))
+
+#' @rdname grad
 setMethod("grad", signature(f = "FunctionObject"), function(f) {
   .FunctionGradCtr(f = f)
 })
 
+#' product overload for FunctionGradObejct
+#'
+#' @param e1 a numeric matrix.
+#' @param e2 a FunctionGradObject created by \code{grad} function.
+#' @return A FunctionGradObject.
+#' @export 
 setMethod("*", signature=c(e1="matrix", e2="FunctionGradObject"),
           function(e1,e2){
         .FunctionGradCtr(f = e2$f, K = e1)    
@@ -38,6 +62,13 @@ setMethod("*", signature=c(e1="matrix", e2="FunctionGradObject"),
 )
 
 ## sum of differential operators
+
+#' plus operator overload for DiffOpObject
+#'
+#' @param e1 a DiffOpObject.
+#' @param e2 a DiffOpObject.
+#' @return A S4 object representing the sum of two differential operators.
+#' @export 
 setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
           function(e1, e2) {
             if (tracemem(e1$f) != tracemem(e2$f)) {
@@ -55,7 +86,13 @@ setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
           }
 )
 
-# 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-diff-op
+#' @export 
 setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
           function(e1, e2) {
             if (tracemem(e1$f) != tracemem(e2$f)) {
@@ -74,7 +111,13 @@ setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"),
               f = e1$f
             )
 })
-## minus (unary) operator for DiffOpObject
+
+#' minus (unary) operator for DiffOpObject
+#'
+#' @param e1 a DiffOpObject.
+#' @return A S4 object representing the sum of two differential operators.
+#' @rdname minus-diff-op
+#' @export 
 setMethod("-", signature(e1 = "DiffOpObject", e2 = "missing"),
           function(e1){
             e1$params[[1]] <- -e1$params[[1]]
@@ -82,7 +125,12 @@ setMethod("-", signature(e1 = "DiffOpObject", e2 = "missing"),
           }
 )
 
-## differential operator product by scalar
+#' 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.
+#' @export 
 setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"),
           function(e1,e2){
             #if (!is.numeric(e1)) stop("bad product")
@@ -96,8 +144,24 @@ setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"),
     Class = "DiffusionOperator",
     contains = "DiffOpObject"
 )
+
 ## laplace() returns a special operator for the case of
 ## isotropic  and stationary diffusion
+
+#' laplace operator
+#'
+#' @param f a FunctionObject.
+#' @return A S4 object representing the laplace operator applied to the function passed as parameter.
+#' @export 
+#' @examples
+#' \dontrun{
+#' library(femR)
+#' data("unit_square")
+#' mesh <- Mesh(unit_square)
+#' Vh <- FunctionSpace(mesh)
+#' f <- Function(Vh)
+#' laplace_f <- laplace(f)
+#' }
 laplace <- function(f) {
     if (!is(f, "FunctionObject")) {
         stop("wrong argument type")
@@ -108,7 +172,24 @@ laplace <- function(f) {
         f = f
     )
 }
+
 ## the general non-isotrpic, non-stationary diffusion operator
+
+#' divergence operator
+#'
+#' @param f a FunctionObject.
+#' @return A S4 object representing the diffusion term of a second order linear differential operator.
+#' @export
+#' @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))
+#' } 
 div <- function(f) {
     if (is(f, "FunctionGradObject")) {
         if (!is.null(f$K)) {
@@ -127,7 +208,27 @@ div <- function(f) {
     Class = "TransportOperator",
     contains = "DiffOpObject"
 )
+
+#' dot product
+#'
+#' @param op1 a numeric vector.
+#' @param op2 a FunctionGradObject.
+#' @return A S4 object representing the advection term of a second order linear differential operator.
+#' @rdname dot-product
+#' @export
+#' @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))
+#' } 
 setGeneric("dot", function(op1, op2) standardGeneric("dot"))
+
+#' @rdname dot-product
 setMethod("dot", signature(op1 = "vector", op2 = "FunctionGradObject"),
           function(op1, op2) {
               .TransportCtr(
@@ -143,6 +244,21 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGradObject"),
     contains = "DiffOpObject"
 )
 
+#' product overload 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.
+#' @export 
+#' @examples
+#' \dontrun{
+#' library(femR)
+#' data("unit_square")
+#' mesh <- Mesh(unit_square)
+#' Vh <- FunctionSpace(mesh)
+#' f <- Function(Vh)
+#' reaction <- 2*f
+#' }
 setMethod("*", signature = c(e1="numeric", e2="FunctionObject"),
           function(e1,e2){
             .ReactionCtr(
@@ -158,6 +274,21 @@ setMethod("*", signature = c(e1="numeric", e2="FunctionObject"),
 )
 
 # overloading stats::dt function
+
+#' dt function overload for FunctionObejct
+#'
+#' @param x a FunctionObject:
+#' @return A S4 object representing the time derivative of a FunctionObject.
+#' @export 
+#' @examples
+#' \dontrun{
+#' library(femR)
+#' data("unit_square")
+#' mesh <- Mesh(unit_square)
+#' Vh <- FunctionSpace(mesh)
+#' f <- Function(Vh)
+#' dt(f)
+#' }
 setMethod("dt", signature = c(x="FunctionObject", df="missing", ncp="missing"), 
           function(x,df,ncp){
   
diff --git a/R/pde.R b/R/pde.R
index 8d54e26..1df0e32 100644
--- a/R/pde.R
+++ b/R/pde.R
@@ -59,8 +59,18 @@
   )
 )
 
-
+#' A PDEs object
+#'
+#' @param L a differential operator.
+#' @param u the forcing term of the PDE.
+#' @param dirichletBC the Dirichlet boundary conditions to be imposed on the boundary of the domain. If not provided, homogeneous Dirichlet boundary conditions will be imposed.
+#' @param initilaCondition the initial condition of a parabolic problem.
+#' @return A S4 object representing a PDE.
+#' @rdname pde
+#' @export 
 setGeneric("Pde", function(L,u,dirichletBC, initialCondition) standardGeneric("Pde"))
+
+#' @rdname pde
 setMethod("Pde", signature=c(L="DiffOpObject", u="ANY", dirichletBC="ANY", 
                              initialCondition="missing"),
           function(L,u,dirichletBC){
@@ -129,6 +139,7 @@ setMethod("Pde", signature=c(L="DiffOpObject", u="ANY", dirichletBC="ANY",
                     is_init=is_init)
           })
 
+#' @rdname pde
 setMethod("Pde", signature=c(L="DiffOpObject", u="ANY", dirichletBC="ANY", 
                              initialCondition="ANY"),
           function(L,u,dirichletBC, initialCondition){
@@ -201,6 +212,7 @@ setMethod("Pde", signature=c(L="DiffOpObject", u="ANY", dirichletBC="ANY",
                     is_init=is_init)
           })
 
+#' @rdname pde
 setMethod("Pde", signature=c(L="DiffOpObject", u="ANY", dirichletBC="missing", 
                              initialCondition="missing"),
           function(L,u,dirichletBC){
diff --git a/R/read_mesh.R b/R/read_mesh.R
index 537ce23..270b398 100644
--- a/R/read_mesh.R
+++ b/R/read_mesh.R
@@ -1,4 +1,10 @@
 # Reads .mesh file (FreeFem++) 
+
+#' reads mesh file from FreeFem++ 
+#'
+#'@param filename path to a .mesh file created by FreeFem++.
+#'@returns a named list to be passed to \code{Mesh}.
+#'@export 
 read_mesh <-function(filename){
   
   # Reading file
diff --git a/src/core b/src/core
index 09c8366..b6b4a5f 160000
--- a/src/core
+++ b/src/core
@@ -1 +1 @@
-Subproject commit 09c8366678cd18a5e7b1b854f1c1fd08703803a2
+Subproject commit b6b4a5f916d92a92981cc975fef1b55bb5b67577