%  back_solve.m    Back substitution with a triangular matrix

function x = back_solve(A,b,display)

%  On entry: A is a square upper triangular matrix,
%            b is a column vector of the same dimension,
%            display is 1 if step-by-step display desired, 0 otherwise.
%  On exit:  x is the solution of Ax=b, assuming A is nonsingular.

n = length(b);
if size(A) ~= [n,n] | size(b) ~= [n,1],
  error('mismatched dimension')
end;

x = zeros(n,1);          % Force x to be a column vector.
x(n) = b(n)/A(n,n);
if display, 
  x, pause,
end;
for i=n-1:-1:1,
  x(i) = (b(i) - A(i,i+1:n)*x(i+1:n))/A(i,i);
  if display,            % Display x after each step if desired.
    x, pause,
  end;
end;