diff --git a/src/chainTDVP.jl b/src/chainTDVP.jl
index 1236c57..bfb9613 100644
--- a/src/chainTDVP.jl
+++ b/src/chainTDVP.jl
@@ -287,7 +287,13 @@ function updatedims(A::Vector, PAs::Vector, PCs::Vector, th, Dlim)
     acc += norm(PAs[N][1:newdims[N],:,1:newdims[N+1]])^2
     return Dims(newdims), effect, acc
 end
-        
+
+"""
+    tdvp1sweep!(dt2, A::Vector, M::Vector, F=nothing; verbose=false, kwargs...)
+
+Propagates the MPS A with the MPO M following the 1-site TDVP method. The sweep is done back and forth with a time step dt2/2. F represents the merged left and right parts of the site being propagated.  
+"""
+
 function tdvp1sweep!(dt2, A::Vector, M::Vector, F=nothing; verbose=false, kwargs...)
     N = length(A)
     dt = dt2/2
diff --git a/src/measure.jl b/src/measure.jl
index 028c0ef..4fe6aaa 100644
--- a/src/measure.jl
+++ b/src/measure.jl
@@ -7,7 +7,24 @@ struct OneSiteObservable <: Observable
     hermitian::Bool
     allsites::Bool
 end
+
+"""
+    OneSiteObservable(name,op,sites)
+
+Computes the local expectation value of the one-site operator `op` on the specified sites. Used to define
+one-site observables that are obs and convobs parameters for the `runsim` function.
+
+"""
 OneSiteObservable(name, op, sites) = OneSiteObservable(name, op, sites, ishermitian(op), false)
+
+"""
+    OneSiteObservable(name,op)
+
+Computes the local expectation value of the one-site operator `op` on the every site. Used to define
+one-site observables that are obs and convobs parameters for the `runsim` function.
+
+"""
+
 OneSiteObservable(name, op) = OneSiteObservable(name, op, nothing, ishermitian(op), true)
 
 struct TwoSiteObservable <: Observable
@@ -33,6 +50,14 @@ end
 CdagCdn(sites::Tuple{Int,Int}) = CdagCdn("CdagCup", sites)
 CdagCdn(i1::Int, i2::Int) = CdagCdn("CdagCdn", (i1,i2))
 
+"""
+    TwoSiteObservable(name,op1,op2,sites1=nothing,sites2=nothing)
+
+Computes the local expectation value of operators `op1` and `op2` where `op1` acts on sites1 and `op2` acts on sites2. Used to define
+several-site observables that are obs and convobs parameters for the `runsim` function.
+
+"""
+
 function TwoSiteObservable(name, op1, op2, sites1=nothing, sites2=nothing)
     return TwoSiteObservable(name, op1, op2, sites1, sites2, sites1==nothing && sites2==nothing)
 end
@@ -68,6 +93,14 @@ reach(A, ob::CdagCup) = max(ob.sites...)
 reach(A, ob::CdagCdn) = max(ob.sites...)
 reach(A, ob::Observable) = 1
 
+"""
+    measurempo(A::Vector, M::Vector)
+
+For a list of tensors `A` representing a right orthonormalized MPS, compute the local expectation
+value of the MPO M on every site.
+
+"""
+
 function measurempo(A::Vector, M::Vector)
     N = length(M)
     N == length(A) || throw(ArgumentError("MPO has $N site while MPS has $(length(A)) sites"))
@@ -77,6 +110,14 @@ function measurempo(A::Vector, M::Vector)
     end
     real(only(F))
 end
+"""
+    measurempo(A::Vector, M::Vector, sites::Tuples{Int,Int})
+
+For a list of tensors `A` representing a right orthonormalized MPS, compute the local expectation
+value of the MPO M on specified sites.
+
+"""
+
 function measurempo(A::Vector, M::Vector, sites::Tuple{Int, Int})
     N = sites[2] - sites[1] + 1
     F = fill!(similar(M[1], (1,1,1)), 1)
@@ -148,9 +189,10 @@ measure(A, O::TwoSiteObservable, ρ::Vector) =
     measure2siteoperator(A, O.op1, O.op2, O.sites1, O.sites2, ρ)
 
 """
-    measure1siteoperator(A, O)
+    measure1siteoperator(A::Vector, O, sites::Vector{Int})
+
 For a list of tensors `A` representing a right orthonormalized MPS, compute the local expectation
-value of a one-site operator O for every site or just one if i is specified.
+value of a one-site operator O for every site or just one if it is specified.
 
 For calculating operators on single sites this will be more efficient if the site is on the left of the mps.
 
@@ -171,6 +213,13 @@ function measure1siteoperator(A::Vector, O, sites::Vector{Int})
     return expval[sites]
 end
 
+"""
+    measure1siteoperator(A::Vector, O, chainsection::Tuple{Int64,Int64})
+For a list of tensors `A` representing a right orthonormalized MPS, compute the local expectation
+value of a one-site operator O for a chainsection.
+
+"""
+
 function measure1siteoperator(A::Vector, O, chainsection::Tuple{Int64,Int64})
     ρ = ones(ComplexF64, 1, 1)
 
@@ -200,6 +249,12 @@ function measure1siteoperator(A::Vector, O, chainsection::Tuple{Int64,Int64})
     return expval
 end
 
+"""
+    measure1siteoperator(A::Vector, O)
+For a list of tensors `A` representing a right orthonormalized MPS, compute the local expectation
+value of a one-site operator O for every site.
+
+"""
 function measure1siteoperator(A::Vector, O)
     N = length(A)
     ρ = ones(ComplexF64, 1, 1)
@@ -214,6 +269,14 @@ function measure1siteoperator(A::Vector, O)
     return expval
 end
 measure1siteoperator(A::Vector, O, ::Nothing) = measure1siteoperator(A, O)
+
+"""
+    measure1siteoperator(A::Vector, O, site::Int)
+For a list of tensors `A` representing a right orthonormalized MPS, compute the local expectation
+value of a one-site operator O for a single site.
+
+"""
+
 function measure1siteoperator(A::Vector, O, site::Int)
     ρ = ones(ComplexF64, 1, 1)
     T = ishermitian(O) ? Float64 : ComplexF64
@@ -404,6 +467,13 @@ function measure2siteoperator_pair(A::Vector, M1, ρ::Vector; conjugate=false)
     return expval + (conjugate ? expval' : transpose(expval)) - dia
 end
 
+"""
+     measure2siteoperator(A::Vector, M1, M2, sites1::Vector{Int}, sites2::Vector{Int})
+
+Caculate expectation of M1*M2 where M1 acts on sites1 and M2 acts on sites2, assumes A is right normalised.
+
+"""
+
 function measure2siteoperator(A::Vector, M1, M2, sites1::Vector{Int}, sites2::Vector{Int})
     if size(M1) == size(M2)
         herm_cis = ishermitian(M1*M2)
diff --git a/src/mpsBasics.jl b/src/mpsBasics.jl
index 18ab21e..dde7aac 100644
--- a/src/mpsBasics.jl
+++ b/src/mpsBasics.jl
@@ -518,14 +518,36 @@ function elementmpo(M, el...)
     c[1,1]
 end
 
+"""
+    apply1siteoperator!(A, O, sites::Vector{Int})
+
+Apply an operator O on the MPS A. O is acting on several sites ::Vector{Int}. The resulting MPS A is the MPS modified by the operator O.
+
+"""
+
 function apply1siteoperator!(A, O, sites::Vector{Int})
     for i in sites
         @tensor R[a,b,s] := O[s,s']*A[i][a,b,s']
         A[i] = R
     end
 end
+
+"""
+    apply1siteoperator!(A, O, sites::Int)
+
+Apply an operator O on the MPS A. O is acting on only one site ::Int. The resulting MPS A is the MPS modified by the operator O.
+
+"""
+
 apply1siteoperator!(A, O, site::Int) = apply1siteoperator!(A, O, [site])
 
+"""
+    applympo!(A, H)
+
+Apply an MPO H on the MPS A. H must have the same number of site than A. The resulting MPS A is the MPS modified by the mpo H.
+
+"""
+
 function applympo!(A, H)
     N = length(H)
     N == length(A) || throw(ArgumentError("MPO has $N site while MPS has $(length(A)) sites"))
diff --git a/src/treeTDVP.jl b/src/treeTDVP.jl
index 265561a..696ade5 100644
--- a/src/treeTDVP.jl
+++ b/src/treeTDVP.jl
@@ -22,6 +22,12 @@ function orthcentersmps(A::TreeNetwork)
     return B
 end
 
+"""
+    physdims(M::TreeNetwork)
+
+Return the physical dimensions of a tree-MPS or tree-MPO `M`.
+"""
+
 function physdims(M::TreeNetwork)
     N = length(M)
     res = Vector{Int}(undef, N)
@@ -262,6 +268,13 @@ end
 
 tdvp1sweep!(dt, A::TreeNetwork, M::TreeNetwork, F=nothing; verbose=false, kwargs...) =
     tdvp1sweep!(dt, A, M, initenvs(A, M, F), findheadnode(A); verbose=verbose, kwargs...)
+
+"""
+    tdvp1sweep!(dt, A::TreeNetwork, M::TreeNetwork, F::Vector, id::Int; verbose=false, kwargs...)
+
+Propagates the tree-MPS A with the tree-MPO M following the 1-site TDVP method. The sweep is done back and forth with a time step dt/2. F represents the merged left and right parts of the site being propagated.  
+"""
+
 function tdvp1sweep!(dt, A::TreeNetwork, M::TreeNetwork, F::Vector, id::Int; verbose=false, kwargs...)
 
     children = A.tree[id].children
@@ -589,6 +602,13 @@ end
 productstatemps(tree::Tree, physdims::Int, Dmax::Int; state=:Vacuum) =
     productstatemps(tree, ntuple(i -> physdims, length(tree)), Dmax; state=state)
 
+"""
+    mpsembed(A::TreeNetwork, Dmax::Int)
+
+Embed tree-MPS `A` in manifold of max bond-dimension `Dmax`
+
+"""
+
 function mpsembed!(A::TreeNetwork, Dmax::Int)
     tree = deepcopy(A.tree)
     pdims = physdims(A)