diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org
index 6a86caf..8a3fd92 100644
--- a/doc/ADI_scheme.org
+++ b/doc/ADI_scheme.org
@@ -1,6 +1,7 @@
 #+TITLE: Adi 2D Scheme
 #+LaTeX_CLASS_OPTIONS: [a4paper,10pt]
 #+LATEX_HEADER: \usepackage{fullpage}
+#+LATEX_HEADER: \usepackage{amsmath}
 #+OPTIONS: toc:nil
 
 * Input
@@ -104,3 +105,96 @@ bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\
 \end{cases}$
 
 [fn:1] $p$ is called =t0_c= inside code
+
+** Finite differences with nodes as cells' centres
+
+*** The explicit FTCS scheme as in PHREEQC
+
+The 1D diffusion equation is:
+
+#+NAME: eqn:1
+\begin{align}
+\frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha \frac{\partial C }{\partial x} \right) \nonumber \\
+   & = \alpha \frac{\partial^2 C}{\partial x^2}
+\end{align}
+
+We discretize it following a Forward Time, Centered Space finite
+difference scheme where the nodes correspond to the centers of a grid
+such as:
+
+[[./grid_pqc.pdf]]
+
+The left boundary is defined on $x=0$ while the center of the first
+cell is in $x=dx/2$, with $dx=L/n$.
+
+We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1
+and assuming constant $\alpha$:
+
+#+NAME: eqn:2
+\begin{equation}\displaystyle
+   \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{i+1}-C^{j+1}_{i}}{\Delta x}-\frac{C^{j+1}_{i}-C^{j+1}_{i-1}}{\Delta x}}{\Delta x}
+\end{equation}
+
+In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on
+the boundaries of each cell (i.e., $(C_{i+1}-C_i)/\Delta x$ on the
+right boundary of the i-th cell and $(C_{i}-C_{i-1})/\Delta x$ on its
+left cell boundary) and then repeat the differentiation to get the
+second derivative of $C$ on the the cell centre $i$.
+
+This discretization works for all internal cells, but not for the
+boundaries. To properly treat them, we need to account for the
+discrepancy in the discretization.
+
+For the first (left) cell, whose center is at $x=dx/2$, we can
+evaluate the left gradient with the left boundary using such distance,
+calling $l$ the numerical value of a constant boundary condition:
+
+#+NAME: eqn:3
+\begin{equation}\displaystyle
+\frac{C_0^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}-
+\frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x}
+\end{equation}
+
+This expression, once developed, yields:
+
+#+NAME: eqn:4
+\begin{align}\displaystyle
+C_0^{j+1} & =  C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\
+          & =  C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right)
+\end{align}
+
+Which leads to the following term applicable to our approach:
+
+#+NAME: eqn:5
+\begin{equation}\displaystyle
+    -C^j_0} = \left[\frac{\alpha \cdot \Delta t}{\Delta x^2}\right] \cdot \left(C^{j+1}_1 + 2l \right) + \left[ -1 - 3 \cdot \frac{\alpha \cdot \Delta t}{\Delta x^2} \right] \cdot C^{j+1}_0
+\end{equation}
+
+In case of constant right boundary, the finite difference of point
+$C_n$ - calling $r$ the right boundary value - is:
+
+#+NAME: eqn:6
+\begin{equation}\displaystyle
+\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}-
+\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x}
+\end{equation}
+
+Which, developed, gives
+#+NAME: eqn:7
+\begin{align}\displaystyle
+C_n^{j+1} & =  C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\
+          & =  C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right)
+\end{align}
+
+If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]]
+becomes zero and we are left with:
+
+
+#+NAME: eqn:8
+\begin{equation}\displaystyle
+C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n)
+\end{equation}
+
+
+
+A similar treatment can be applied to the BTCS implicit scheme.
diff --git a/doc/grid_pqc.pdf b/doc/grid_pqc.pdf
new file mode 100644
index 0000000..b970d18
Binary files /dev/null and b/doc/grid_pqc.pdf differ