function x = Solve(A,pivots,b)
%  Solve solves the linear system A*x = b using the decomposition
%  obtained from Factor.  Do not use Solve if a zero pivot was
%  detected in Factor, i.e., flag was positive.
%
%  Input variables:
%     A           = output of Factor.  Triangular decomposition
%                   of the coefficient matrix.
%     pivots      = output of Factor. Record of row interchanges.
%     b           = right hand side vector b.
%  Output variables:
%     x           = solution vector x.

b = b(:);
neq = length(b);

if neq > 1
  %  Forward elimination.
  for k = 1:neq-1
    m = pivots(k);
    b([m k]) = b([k m]);
    b(k+1:neq) = b(k+1:neq) + A(k+1:neq,k)*b(k);
  end

  %  Back substitution.
  b(neq) = b(neq) / A(neq,neq);
  for i = neq-1:-1:1
    b(i) = (b(i) - A(i,i+1:neq)*b(i+1:neq)) / A(i,i);
  end
else
  b(1) = b(1)/A(1,1);
end

x = b;
