diff --git a/.gitignore b/.gitignore
index cbb38e68280b1b395ec831e9bd3c0d1c158a258b..93fca5e5d252a76785b7d8a1fef3b9d5d5d5939a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,3 +14,5 @@ GTAGS
 /.clangd/
 /.cache/
 /doc/lisp/elpa/
+/doc/*.png
+/doc/*.gnu
diff --git a/doc/userdoc.org b/doc/userdoc.org
index 15175e2f8180c9c87a35ecf115a2c4153c47b03d..8764d39d9191f2b752a274fa331759117d7b37e1 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -9,7 +9,7 @@
 #+OPTIONS: tex:t
 
 #+LANGUAGE: en
-
+#+ATTR_LATEX: :width 4cm
 #+HTML_HEAD_EXTRA: <style> pre.src-pugs:before { content: 'Pugs'; } </style>
 
 #+LATEX_CLASS_OPTIONS: [10pt]
@@ -28,9 +28,96 @@
 
 * Introduction
 
-Reading this document one should get how to use ~pugs~.
+~pugs~[fn:pugs-def] is a general purpose solver collection designed to
+approximate solutions of partial differential equations. It is mainly
+(but not only) designed to deal with hyperbolic problems using
+finite-volume methods.
+
+~pugs~ is a parallel software that uses ~MPI~[fn:MPI-def] and multi-threading
+for parallelism. Multi-threading is achieved through an encapsulation
+of some [[https:github.com/kokkos/kokkos][Kokkos]] mechanisms.
+
+The philosophy of ~pugs~ is to provide "simple" numerical tools that are
+assembled together through a high-level language (a DSL[fn:DSL-def]
+close to the mathematics) to build more complex solvers. This approach
+is inspired by the success of [[http://freefem.org][FreeFEM]], which use a similar approach.
+
+Before detailing the leading concepts and choices we have made to
+develop ~pugs~, we give a simple example that should illustrate the
+capabilities of the code.
+
+** An example
+For instance, the following code builds a uniform Cartesian grid of
+$]-1,1[^2$ made of $20\times20$ cells and transforms it by displacing the
+mesh nodes according to a user defined vector field $T: \mathbb{R}^2
+\to \mathbb{R}^2$.
+#+NAME: introduction-example
+#+BEGIN_SRC pugs :exports both :results output
+  import mesh;
+  import math;
+  import writer;
+
+  let pi:R, pi = acos(-1);
+  let theta:R^2 -> R, x -> 0.5*pi*(x[0]+1)*(x[0]-1)*(x[1]+1)*(x[1]-1);
+  let M:R^2 -> R^2x2, x -> (cos(theta(x)), -sin(theta(x)),
+                            sin(theta(x)),  cos(theta(x)));
+  let T: R^2 -> R^2, x -> x + M(x)*x;
+  let m:mesh, m = cartesian2dMesh((-1,-1), (1,1), (20,20));
 
-* Basics
+  m = transform(m, T);
+
+  write_mesh(gnuplot_writer("transformed", 0), m);
+#+END_SRC
+The example is quite easy to read:
+- first, some *modules* are loaded: the ~mesh~ module, which contains
+  some mesh manipulation functions. The ~math~ module which provides a
+  set of classical mathematical functions ($\sin$, $\cos$,
+  ...). The ~writer~ module is used to save meshes or discrete functions
+  to files using various formats.
+- The second block of data defines variables of different kind
+  - ~pi~ is a real value that is initialized by an approximation of $\pi$.
+  - ~theta~ is the real value function $\theta$ defined by
+        $$\theta: \mathbb{R}^2 \to \mathbb{R},\quad\mathbf{x} \mapsto \frac{\pi}{2} (x_0+1)(x_0-1)(x_1+1)(x_1-1) $$
+  - ~M~ is the $\mathbb{R}^{2\times2}$  matrix field $M$ defined by
+        $$M: \mathbb{R}^2 \to \mathbb{R}^{2\times2},\quad\mathbf{x} \mapsto \begin{pmatrix}
+     \cos(\theta(\mathbf{x})) & -\sin(\theta(\mathbf{x}))  \\
+     \sin(\theta(\mathbf{x}))  & \cos(\theta(\mathbf{x}))
+     \end{pmatrix}$$
+  - ~T~ is the vector field $T$ defined by
+    $$
+    T: \mathbb{R}^2 \to \mathbb{R}^2, \quad\mathbf{x} \mapsto (I+M(\mathbf{x}))\mathbf{x}
+    $$
+  - Finally ~m~ is defined as the uniform Cartesian mesh grid of
+    $]-1,1[^2$. The last argument: ~(20,20)~ sets the number of cells in
+    each direction: $20$.
+- The third block is the application of the transformation $T$ to the
+  mesh. Observe that if the resulting mesh is stored in the same
+  variable ~m~, the old one was not modified in the process. It is
+  important to already have in mind that the ~pugs~ language *does not
+  allow* the modifications of values of *non-basic* types. This is
+  discussed in the section [[basics]].
+- Finally, the last block consists in saving the obtained mesh in a
+  ~gnuplot~ file. The result is shown on Figure [[fig:intro-example]].
+
+#+NAME: intro-transform-mesh-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/intro-transform-mesh.png")
+  reset
+  unset grid
+  unset border
+  unset key
+  unset xtics
+  unset ytics
+  set size square
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/transformed.0000.gnu)' w l
+#+END_SRC
+
+#+CAPTION: Obtained transformed mesh
+#+NAME: fig:intro-example
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: intro-transform-mesh-img
+
+*  Basics<<basics>>
 
 ** Hello world!
 
@@ -41,7 +128,7 @@ Reading this document one should get how to use ~pugs~.
 
 *** Pugs
 In this simple example, one defines the function $f:\mathbb{R} \to \mathbb{R}, x \mapsto 2\sin(x)$
-#+name: hello-world
+#+NAME: hello-world
 #+BEGIN_SRC pugs :exports both :results output
   import math;
   let f: R -> R, x -> 2*sin(x);
@@ -49,9 +136,13 @@ In this simple example, one defines the function $f:\mathbb{R} \to \mathbb{R}, x
   cout << "f(12) = " <<  f(12) << "\n";
 #+END_SRC
 Then one prints ~Hello world!~ and the evaluation of $f$ at position $12$.
-#+results: hello-world
+#+RESULTS: hello-world
 * The end
 
 #+BEGIN_SRC python :exports both :results output
   print ("hello world!")
 #+END_SRC
+
+[fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers
+[fn:MPI-def] ~MPI~: Message Passing Interface
+[fn:DSL-def] ~DSL~: Domain Specific Language