diff --git a/README.md b/README.md
index cab2b2f..1ad3707 100644
--- a/README.md
+++ b/README.md
@@ -8,8 +8,10 @@
 
 This package is inspired by [SYCOMORE](https://sycomore.readthedocs.io/) with the intention to use Automatic differentiation (with ForwardDiff.jl)
 
+Take a look at the [documentation](https://atrotier.github.io/EPGsim.jl/dev/).
+
 ## TO DO list :
-- Delete states with amplitude < defined threshold
--  discrete / discrete 3D
+- Delete states with amplitude < defined threshold or number of states (implemented in epgDephasing ??)
+- discrete / discrete 3D
 - Diffusion
-- Magnetization transfer
\ No newline at end of file
+- Magnetization Transfer
\ No newline at end of file
diff --git a/docs/src/AD.md b/docs/src/AD.md
index 4ab4069..31c9cda 100644
--- a/docs/src/AD.md
+++ b/docs/src/AD.md
@@ -111,6 +111,36 @@ lines!(ax,TE_vec,df)
 f
 ```
 
+## Differentiation along multiple variables
+If we want to obtain the derivation along T1 and T2 we need to change the EPG_MESE function. The function should take as input a vector containing T2 and T1 (here noted T2/T1) :
+```@example AD
+function MESE_EPG2(T2T1,TE,ETL,delta)
+  T2,T1 = T2T1
+  T = complex(eltype(T2))
+  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
+  echo_vec = Vector{Complex{eltype(T2)}}()
+
+  E = epgRotation(E,pi/2*delta, pi/2)
+  ##loop over refocusing-pulses
+  for i = 1:ETL
+    E = epgDephasing(E,1)
+    E = epgRelaxation(E,TE,T1,T2)
+    E  = epgRotation(E,pi*delta,0.0)
+    E  = epgDephasing(E,1)
+    push!(echo_vec,E.Fp[1])
+  end
+
+  return abs.(echo_vec)
+end
+
+j2 = ForwardDiff.jacobian(x -> MESE_EPG2(x,TE,ETL,deltaB1),[T2,T1])
+```
+Here we can see that the second column corresponding to T1 is equal to 0 which is expected for a MESE sequence and the derivative along T2 gives the same results :
+
+```@example AD
+j2[:,1] ≈ vec(j)
+```
+
 ## Reproducibility
 
 This page was generated with the following version of Julia:
diff --git a/docs/src/index.md b/docs/src/index.md
index e72077c..53e4c1a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -14,7 +14,7 @@ The principal aspect of this package was to make it compatible with Automatic Di
 !!! note
     EPGsim.jl is work in progress and in some parts not entirely optimized. The interface is susceptible to change between version
 
-## physics concept
+## EPG concept 
 Introduction to the physics concepts behing EPG as well as their usage can be found on the rad229 youtube channels by Brian Hargreaves and Daniel Ennis :
 - [Lecture-04A: Definition of the Extended Phase Graph Basis](https://www.youtube.com/watch?v=bskhnaoJVNY)
 - [Lecture-04B: Sequence Operations in the Extended Phase Graph Domain](https://www.youtube.com/watch?v=kToL-9ZTzCs)
diff --git a/docs/src/regular.md b/docs/src/regular.md
index 8ef7076..311e08b 100644
--- a/docs/src/regular.md
+++ b/docs/src/regular.md
@@ -54,6 +54,10 @@ E = epgDephasing(E,1)
 E = epgRotation(E,deg2rad(60),deg2rad(117))
 ```
 
+!!! note
+    Currently, all the EPGstates are stored and used for calculation. 
+    The states equal or really close to zero are not deleted
+
 # Accessing states
 States can seen directly as a vector :
 ```@example Regular
diff --git a/src/RegularEPG.jl b/src/RegularEPG.jl
index 86a03f6..e5d4bb5 100644
--- a/src/RegularEPG.jl
+++ b/src/RegularEPG.jl
@@ -1,6 +1,23 @@
 export epgDephasing,epgRelaxation,epgRotation,rfRotation
 export EPGStates, getStates
 
+"""
+    EPGStates{T <: Real} 
+
+Stores the EPG states in 3 vectors Fp,Fn and Z.
+
+# Constructors :
+    EPGStates(Fp::Vector{Complex{S}},Fn::Vector{Complex{S}},Z::Vector{Complex{S}}) where {S <: Real}
+    EPGStates(Fp::T=0,Fn::T=0,Z::T=1) where T <: Number
+
+# Fields
+- `Fp::Vector{Complex{T}}`
+- `Fn::Vector{Complex{T}}`
+- `Z::Vector{Complex{T}}`
+
+# Related functions
+- `getStates(E::EPGStates)` : extract EPG states as matrix 3xN
+"""
 mutable struct EPGStates{T <: Real} 
   Fp::Vector{Complex{T}}
   Fn::Vector{Complex{T}}
@@ -18,7 +35,11 @@ mutable struct EPGStates{T <: Real}
   end
 end
 
+"""
+    getStates(E::EPGStates)
 
+Extract EPG states as matrix 3xN
+"""
 function EPGStates(Fp::T=0,Fn::T=0,Z::T=1) where T <: Number
   T2 = ComplexF64
 
@@ -33,8 +54,9 @@ end
     epgDephasing(E::EPGStates, n=1) where T
   
 shifts the transverse dephasing states `F` corresponding to n dephasing-cycles.
+n can be any integer
 """
-function epgDephasing(E::EPGStates, n=1)
+function epgDephasing(E::EPGStates, n::Int=1)
   
   if(abs(n)>1)
     for i in 1:abs(n)
@@ -88,9 +110,13 @@ function epgRelaxation(E::EPGStates,t,T1, T2)
 end
 
 """
-    rfRotation_AT(alpha, phi=0.)
+    rfRotation(alpha, phi=0.)
 
 returns the rotation matrix for a pulse with flip angle `alpha` and phase `phi`.
+
+  # Arguments
+* `alpha`  - flip angle (radian)
+* `phi=0.` - phase of the flip angle (radian)
 """
 function rfRotation(alpha, phi=0.)
   R = [ cos(alpha/2.)^2   exp(2*im*phi)*sin(alpha/2.)^2   -im*exp(im*phi)*sin(alpha);