Suppose that we are solving a system of *n* linear equations
**L x** = **b** where **L** = (l_{ij}) is a lower-triangular
matrix with no zero entries on the diagonal. The steps in
solving for **x** = (*x*_{i}) are:

For *i* = 1, 2, ..., *n*, in that order, let:

This is called forward substitution because we are starting
with the first unknown *x*_{1} and, having solved for
it, using it to solve for the next unknown *x*_{2}, and so on.
In Matlab, this may be done by the (poorly implemented) code:

n = length( b );
x = zeros( n, 1 );
for i=1:n
x(i) = b(i);
for j=1:(i - 1)
x(i) = x(i) - L(i, j)*x(j);
end
x(i) = x(i)/L(i, i);
end

Note that if we initialize the vector **x** to be the zero vector,
then the sum is simply the dot product of the *i*th row of A
dotted with *x*. Thus, we may code this algorithm in Matlab
as:

n = length( b );
x = zeros( n, 1 );
for i=1:n
x(i) = ( b(i) - L(i, :)*x )/L(i, i);
end

Note: this is **significantly faster** than the previous implementation
(and is also easier to remember).

In the special case where the diagonal entries of **L** are all
zero, the above Matlab code simplifies to:

n = length( b );
x = zeros( n, 1 );
for i=1:n
x(i) = b(i) - L(i, :)*x;
end