the channel at next time step, : initial condition for the previous channel to the one analyzed at present time step, : initial condition for the next channel to the one analyzed at present time step. For the last cell of the channels (z), the most convenient condition becomes an adiabatic one, so This condition is satisfied if and only if (regardless of a null value) Let us solve this problem (in a matrix form) for the case of 3 channels and 5 nodes (including the initial boundary condition). We express this as a linear system problem: where Now we must realize that AA and BB should be arrays made of four different subarrays (remember that only three channels are considered for this example, but it covers the main part discussed above): where the elements mentioned above correspond to the next arrays, and an additional 4×4 full of zeros. Please note that the sizes of AA and BB are 12×12: The d vector here is used to hold the boundary conditions. In this example it is a 12×1 vector: To find the concentration at any time, one must iterate the following equation: Example: 2D diffusion When extending into two dimensions on a uniform Cartesian grid, the derivation is similar and the results may lead to a system of band-diagonal equations rather than tridiagonal ones. The two-dimensional heat equation can be solved with the Crank–Nicolson discretization of assuming that a square grid is used, so that . This equation can be simplified somewhat by rearranging terms and using the CFL number For the Crank–Nicolson numerical scheme, a low CFL number is not required for stability, however, it is required for numerical accuracy. We can now write the scheme as Solving such a linear system is not practical due to extremely high time complexity of solving a linear system by the means of Gaussian elimination or even Strassen algorithm. Hence an alternating-direction implicit method can be implemented to solve the numerical PDE, whereby one dimension is treated implicitly, and other dimension explicitly for half of the assigned time step and conversely for the remainder half of the time step. The benefit of this strategy is that the implicit solver only requires a tridiagonal matrix algorithm to be solved. The difference between the true Crank–Nicolson solution and ADI approximated solution has an order of accuracy of and hence can be ignored with a