## 2 Development of the Procedure  
The easiest way to discretize a partial differential equation problem (like we have on Equation (1)) consist on reducing the continuous domain into an equispaced finite number of points \((x_{i},t_{j})\) located on the nodes of a uniform rectangular mesh. The partial derivatives of \(u\) can be estimated now by the following difference quotients  
\[\frac{\partial u(x_{i},t_{j})}{\partial x}\approx\frac{u_{i+1,j}-u_{i,j}}{h}, \qquad\frac{\partial u(x_{i},t_{j})}{\partial t}\approx\frac{u_{i,j+1}-u_{i, j}}{k} \tag{7}\]  
which are called forward differences and have a truncation error of order 1. The partial derivatives of \(u\) can also be obtained by  
\[\frac{\partial u(x_{i},t_{j})}{\partial x}\approx\frac{u_{i,j}-u_{i-1,j}}{h}, \qquad\frac{\partial u(x_{i},t_{j})}{\partial t}\approx\frac{u_{i,j}-u_{i,j-1}}{k} \tag{8}\]  
which are backward differences and have also a truncation error of order 1. Finally the partial derivatives of \(u\) can also be obtained by the following  
\[\frac{\partial u(x_{i},t_{j})}{\partial x}\approx\frac{u_{i+1,j}-u_{i-1,j}}{2 h},\qquad\frac{\partial u(x_{i},t_{j})}{\partial t}\approx\frac{u_{i,j+1}-u_{i,j-1}}{2k} \tag{9}\]  
\[\frac{\partial^{2}u(x_{i},t_{j})}{\partial x^{2}}\approx\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}},\qquad\frac{\partial^{2}u(x_{i},t_{j})}{\partial t^{2}} \approx\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{k^{2}} \tag{10}\]  
which are called central differences and have a truncation error of order 2. In these expressions, \(u_{i,j}\) denotes the approximated value of the unknown \(u\) function on \((x_{i},t_{j})\), \(h=\frac{b-a}{n}\) and \(k=\frac{t_{max}-t_{min}}{m}\) are respectively the space and time step on the mesh and \(n\) and \(m\) are respectively the space and time amount of subintervals to be considered. Then, we consider \(x_{i+1}=x_{i}+h\) and \(t_{j+1}=t_{j}+k\), where \(i=0,1,...,n-1\) and \(j=0,1,...,m-1\).  
Let us apply Crank-Nicolson difference scheme, which consists of discretizing Equation (1) in two consecutive time instants of \(j\) and average them. In the first instant \(t_{j}\) we approximate \(\frac{\partial u}{\partial t}\) by the forward difference described on Equation (7) and in the second instant \(t_{j+1}\) we do it by the backward one. At both instants we approximate \(\frac{\partial u}{\partial x}\) and \(\frac{\partial^{2}u}{\partial x^{2}}\) by using the central differences (described by Equations (9) and (10) respectively). By means of this procedure we get, for the first instant,  
\[\frac{u_{i,j+1}-u_{i,j}}{k}-\varepsilon\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{ 2}}+u_{i,j}\frac{u_{i+1,j}-u_{i-1,j}}{2h}=0 \tag{11}\]  
and for the second one,  
\[\frac{u_{i,j+1}-u_{i,j}}{k}-\varepsilon\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1 }}{h^{2}}+u_{i,j+1}\frac{u_{i+1,j+1}-u_{i-1,j+1}}{2h}=0 \tag{12}\]  
After averaging both expressions, it results  
\[Bu_{i,j+1}u_{i+1,j+1}-Bu_{i,j+1}u_{i-1,j+1}-Cu_{i+1,j+1}+2 Cu_{i,j+ 1}-Cu_{i-1,j+1}+Au_{i,j+1}= \tag{13}\] \[Cu_{i+1,j}+(A-2C)u_{i,j}+Cu_{i-1,j}+Bu_{i,j}u_{i+1,j}-Bu_{i,j}u_{ i-1,j}\]  
where \(A=\frac{1}{k}\), \(B=\frac{1}{4h}\), \(C=\frac{\varepsilon}{2h^{2}}\) and the range of indices is \(i=1,...,n-1\) and \(j=0,1,...,m-1\). This expression (which has a truncation error of order 2) results on a \(m\) nonlinear systems of \(n-1\) unknowns and \(n-1\) equations. Each one of these systems comes after considering that we already know the value of \(u_{0,j}\) and \(u_{n,j}\), \(\forall j\) (from the boundary conditions) and \(u_{i,0}\), \(\forall i\) from the initial condition. Let us notice that the values corresponding to the \(j\) subindex have been calculated on the previous systems and the only unknowns will be those with the \(j+1\) subindex.