We can rewrite our original differential equation in terms of a discretized equation with approximations to the derivatives as -\frac{u_{i+1} -2u_i +u_{i-i}}{h^2}=f(x_i,u(x_i)), with i=1,2,\dots, n . We need to add to this system the two boundary conditions u(a) =u_0 and u(b) = u_{n+1} . If we define a matrix \mathbf{A} = \frac{1}{h^2}\begin{bmatrix} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & \dots & \dots &\dots &\dots & \dots \\ & & &-1 &2& -1 \\ & & & &-1 & 2 \\ \end{bmatrix} and the corresponding vectors \mathbf{u} = (u_1, u_2, \dots,u_n)^T and \mathbf{f}(\mathbf{u}) = f(x_1,x_2,\dots, x_n,u_1, u_2, \dots,u_n)^T we can rewrite the differential equation including the boundary conditions as a system of linear equations with a large number of unknowns \mathbf{A}\mathbf{u} = \mathbf{f}(\mathbf{u}).