diff --git a/doc/src/SupportVMachines/svm.do.txt b/doc/src/SupportVMachines/svm.do.txt
index f16bcfadd..a645290ed 100644
--- a/doc/src/SupportVMachines/svm.do.txt
+++ b/doc/src/SupportVMachines/svm.do.txt
@@ -21,10 +21,13 @@ definition of a _margin_ which separates classes (in case of
 classification problems) of variables. It is also used for regression
 problems.
 
-With SVMs we distinguish between hard margin and soft margins. The latter introduces a so-called softening parameter to be discussed below.
-We distinguish also between linear and non-linear approaches. The latter are the most frequent ones since it is rather unlikely that we can separate classes easily by say straight lines. 
-
+With SVMs we distinguish between hard margin and soft margins. The
+latter introduces a so-called softening parameter to be discussed
+below.  We distinguish also between linear and non-linear
+approaches. The latter are the most frequent ones since it is rather
+unlikely that we can separate classes easily by say straight lines.
 
+_Note: several figures are missing. They will be added shortly. To run the codes, use the jupyter notebook_
 
 
 !split
@@ -553,17 +556,17 @@ we need to introduce for example a polynomial transformation to a two-dimensiona
 !split
 ===== The equations =====
 
-Suppose we define a polynomial transformation of degree two (we continue to live in a plane with $x_1$ and $x_2$ as variables)
+Suppose we define a polynomial transformation of degree two only (we continue to live in a plane with $x_i$ and $y_i$ as variables)
 !bt
 \[
-z = \phi(x) =\left(1, x_1, x_2, x_1^2, x_2^2, x_1x_2).
+z = \phi(x_i) =\left(x_i^2, y_i^2, \sqrt{2}x_iy_i\right).
 \]
 !et
 
 With our new basis, the equations we solved earlier are basically the same, that is we have now (without the slack option for simplicity)
 !bt
 \[
-{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{z}_i^T\bm{Z}_j,
+{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{x}_i^T\bm{z}_j,
 \]
 !et
 subject to the constraints $\lambda_i\geq 0$, $\sum_i\lambda_iy_i=0$, and for the support vectors
@@ -572,30 +575,95 @@ subject to the constraints $\lambda_i\geq 0$, $\sum_i\lambda_iy_i=0$, and for th
 y_i(\bm{w}^T\bm{z}_i+b)= 1 \hspace{0.1cm}\forall i,
 \]
 !et
-from which we also find $b$. 
+from which we also find $b$.
+To compute $\bm{z}_i^T\bm{Z}_j$ we define the kerne $K(\bm{x}_i,\bm{x}_j)$ as
+!bt
+\[
+K(\bm{x}_i,\bm{x}_j)=\bm{z}_i^T\bm{Z}_j= \phi(\bm{x}_i)^T\phi(\bm{x}_j).
+\]
+!et
+For the above example, the kernel reads
+!bt
+\[
+K(\bm{x}_i,\bm{x}_j)=[x_i^2, y_i^2, \sqrt{2}x_iy_i]^T\begin{bmatrix} x_j^2 \\ y_j^2 \\ \sqrt{2}x_jy_j \end{bmatrix}=x_i^2x_j^2+2x_ixJy_iy_j+y_i^2y_j^2.
+\]
+!et
+We note that this is nothing but the dot product of the two original vectors $(\bm{x}_i^T\bm{x}_j)^2$. Instead of thus computing the product in the Lagrangian of $\bm{x}_i^T\bm{z}_j$ we simply compute the dot product $(\bm{x}_i^T\bm{x}_j)^2$. This leads to the so-called kernel trick and the result leads to the same as if we went through the trouble of performing the transformation $(\phi(\bm{x}_i)^T\phi(\bm{x}_j)$ during the SVM calculations. 
+
 
 !split
-===== Different kernels =====
+===== The problem to solve =====
+Using our definition of the kernel We can rewrite again the Lagrangian
+!bt
+\[
+{\cal L}=\sum_i\lambda_i-\frac{1}{2}\sum_{ij}^n\lambda_i\lambda_jy_iy_j\bm{x}_i^T\bm{z}_j,
+\]
+!et
+subject to the constraints $\lambda_i\geq 0$, $\sum_i\lambda_iy_i=0$ in terms of a convex optimization problem
+!bt
+\[
+\frac{1}{2} \bm{\lambda}^T\begin{bmatrix} y_1y_1K(\bm{x}_1,\bm{x}_1) & y_1y_2K(\bm{x}_1,\bm{x}_2) & \dots & \dots & y_1y_nK(\bm{x}_1,\bm{x}_n) \\
+y_2y_1\bm{x}_2^T\bm{x}_1 & y_2y_2\bm{x}_2^T\bm{x}_2 & \dots & \dots & y_1y_n\bm{x}_2^T\bm{x}_n \\
+\dots & \dots & \dots & \dots & \dots \\
+\dots & \dots & \dots & \dots & \dots \\
+y_ny_1K(\bm{x}_n,\bm{x}_1) & y_ny_2K(\bm{x}_n\bm{x}_2) & \dots & \dots & y_ny_nK(\bm{x}_n,\bm{x}_n) \\
+\end{bmatrix}\bm{\lambda}-\mathbb{I}\bm{\lambda}, 
+\]
+!et
+subject to $\bm{y}^T\bm{\lambda}=0$. Here we defined the vectors $\bm{\lambda} =[\lambda_1,\lambda_2,\dots,\lambda_n]$ and 
+$\bm{y}=[y_1,y_2,\dots,y_n]$. 
+If we add the slack constants this leads to the additional constraint $0\leq \lambda_i \leq C$.
+
+We can rewrite this (see the solutions below) in terms of a convex optimization problem of the type
+!bt
+\begin{align*}
+    &\mathrm{min}_{\lambda}\hspace{0.2cm} \frac{1}{2}\bm{\lambda}^T\bm{P}\bm{\lambda}+\bm{q}^T\bm{\lambda},\\ \nonumber
+    &\mathrm{subject to} \hspace{0.2cm} \bm{G}\bm{\lambda} \preceq \bm{h} \hspace{0.2cm} \bm{A}\bm{\lambda}=f.
+\end{align*}
+!et
+Below we discuss how to solve these equations. Here we note that the matrix $\bm{P} has matrix elements $p_{ij}=y_iy_jK(\bm{x}_i,\bm{x}_j)$.
+Given a kernel $K$ and the targets $y_i$ this matrix is easy to set up. The constraint $\bm{y}^T\bm{\lambda}=0$ leads to $f=0$ and $\bm{A}=\bm{y}$. How to set up the matrix $\bm{G}$ is discussed later. Here note that the inequalities $0\leq \lambda_i \leq C$ can be split up into
+$0\leq \lambda_i$ and $\lambda_i \leq C$. These two inequalities define then the matrix $\bm{G}$ and the vector $\bm{h}$.
+
 
 !split
-===== Quadratic coefficient matrix =====
+===== Different kernels and Mercer's theorem =====
+
+There are several popular kernels being used. These are
+o Linear: $K(\bm{x},\bm{y})=\bm{x}^T\bm{y}$,
+o Polynomial: $K(\bm{x},\bm{y})=(\bm{x}^T\bm{y}+\gamma)^d$,
+o Gaussian Radial Basis Function: $K(\bm{x},\bm{y})=\exp{\left(-\gamma\vert\vert\bm{x}-\bm{y}\vert\vert^2\right)}$,
+o Tanh: $K(\bm{x},\bm{y})=\tanh{(\bm{x}^T\bm{y}+\gamma)}$,
+and many other ones.
+
+An important theorem for us is "Mercer's theorem":https://en.wikipedia.org/wiki/Mercer%27s_theorem". 
+The theorem states that if a kernel function $K$ is symmetric, continuous and leads to a positive definite matrix $\bm{P}$ then
+there exists a function $\phi$ that maps $\bm{x}_i$ and $\bm{x}_j$  into another space
+(possibly with much higher dimensions) such that
+!bt
+\[
+K(\bm{x}_i,\bm{x}_j)=\phi(\bm{x}_i)^T\phi(\bm{x}_j).
+\]
+!et
+So you can use $K$ as a kernel since you know $\phi$ exists, even if
+you don’t know what $\phi$ is. 
+Note that some frequently used kernels (such as the Sigmoid kernel) don’t respect all of Mercer’s conditions, yet they generally work
+well in practice.
 
-!split 
-===== Mercer's theorem =====
 
 !split
 =====  Mathematical optimization of convex functions =====
 
-A mathematical optimization problem, or just optimization problem, has the form
+A mathematical (quadratic) optimization problem, or just optimization problem, has the form
 !bt
-\[
-\mathrm{minimize}\hspace{0.1cm} f(x),
-\]
+\begin{align*}
+    &\mathrm{min}_{\lambda}\hspace{0.2cm} \frac{1}{2}\bm{\lambda}^T\bm{P}\bm{\lambda}+\bm{q}^T\bm{\lambda},\\ \nonumber
+    &\mathrm{subject to} \hspace{0.2cm} \bm{G}\bm{\lambda} \preceq \bm{h} \hspace{0.2cm} \bm{A}\bm{\lambda}=f.
+\end{align*}
 !et
-subject to some constraints $g(\lambda_i) \leq  b_i$ for say a selected set $i=1,2,\dots, n$.
+subject to some constraints for say a selected set $i=1,2,\dots, n$.
 In our case we are optimizing with respect to the Lagrangian multipliers $\lambda_i$, and the
 vector $\bm{\lambda}=[\lambda_1, \lambda_2,\dots, \lambda_n]$ is the optimization variable we are dealing with.
-and $f(x)$ is our objective function while $g(\lambda_i) \leq  b_i$ represents our constraint function.
 
 In our case we are particularly interested in a class of optimization problems called convex optmization problems. 
 In our disussion on gradient descent methods we discussed at length the definition of a convex function. 
@@ -615,19 +683,37 @@ problem using say standard gradient methods or conjugate gradient
 methods. However, these methods tend to exhibit a rather slow
 converge. So, welcome to the promised land of quadratic programming.
 
-The functions we need are contained in the quadratic programming package _CVXOPT_ and we need to import it
+The functions we need are contained in the quadratic programming package _CVXOPT_ and we need to import it together with _numpy_ as
+
 !bc pycod
 import numpy
 import cvxopt
 !ec
 
-Let us first set up the standard form the of quadratic programming (QP) equations by defining the problem as
+This will make our life much easier. You don't need t write your own optimizer.
+
+
+!split
+===== A simplex example =====
+
+Let us show how to perform the optmization using a simple case. Assume we want to optimize the following problem
+!bt
+\begin{align*}
+    &\mathrm{min}_{\lambda}\hspace{0.2cm} \frac{1}{2}x^2+5x+3y \\ \nonumber
+    &\mathrm{subject to} \\ \nonumber
+    &x, y \geq 0 \\ \nonumber
+    &x+3y  \geq 15 \\ \nonumber
+    &2x+5y  \leq  100 \\ \nonumber
+    &3x+4y  \leq  80.  \\ \nonumber
+\end{align*}
+!et
+The minimization problem can be rewritten in terms of vectors and matrices as (with $x$ and $y$ being the unknowns)
 !bt
 \[
-\mathrm{min} 
+frac{1}{2}\begin{bmatrix} x\\ y \end{bmatrix}^T   \begin{bmatrix} 1 & 0\\ 0 & 0  \end{bmatrix}  \begin{bmatrix} x \\ y \end{bmatrix}  + \begin{bmatrix}3\\ 4  \end{bmatrix}^T \begin{bmatrix}x \\ y  \end{bmatrix}.  
 \]
 !et
-
+Similarly, we an now 
 subject to Gx  u. Note that x itself is not provided to the solver, since it is an internal
 variable being optimized over. In particular, this means that the solver has no explicit knowledge
 of x itself; everything is implicity defined by the supplied parameters. It is essential
@@ -666,3 +752,8 @@ sol = solvers.qp(P,q,G,h)
 sol[’x’] # [7.13e-07, 5.00e+00]
 sol[’primal objective’]
 !ec
+
+
+!split
+===== Multiclass problems and regression with SVMs =====
+This material will be added later.