function [x] = nluSolve(A,b)
% nluSolve solve linear system using LU factorization
% Syntax: x = nluSolve(A,b)
% A: n by n matrix returned in LU form by nluFactor
% b: right side vector
% x: solution vector

[m,n] = size(A);
if m ~= n
   error('Matrix is not square');
end

% Do forward substitution to solve Ly = b

y = zeros(n,1);
y(1) = b(1);
for i = 2 : n
   s = b(i);
   for j = 1 : i-1
      s = s - A(i,j)*y(j);
   end
   y(i) = s;
end

% Do backward substitution to solve Ux = y

x = zeros(n,1);
x(n) = y(n)/A(n,n);
for i = n-1: -1 : 1
   s = y(i);
   for j = i+1 : n
      s = s - A(i,j)*x(j);
   end
   x(i) = s / A(i,i);
end

end % nluSolve