function x = substitute ( l,u,b )
%
%  x = substitute(l,u,b) solves l*y = b and u*x = y,
%    where l is lower and u upper triangular.
%
%  example: a = randn(3,3); b = randn(3,1);
%           [l,u] = simple_lu(a);
%           x = substitute(l,u,b);
%           b - a*x
%
n = size(l,1);
y = b;                 % solve l*y = b
for i=2:n
  for j=1:i-1
    y(i) = y(i) - l(i,j)*y(j);
  end;
end;
x = y;                 % solve u*x = y
for i=n:-1:1
  for j=i+1:n
    x(i) = x(i) - u(i,j)*x(j);
  end;
  x(i) = x(i)/u(i,i);
end;