function x = lu_solv(LU,b)
% Solves A*x = L*U*x = b, where LU is obtained by LU decomposition of A and
% contains both L and U as [L\U]. b must be a column vector.

n = length(b);

for i = 2:n
    b(i) = b(i) - LU(i,1:i-1)*b(1:i-1);
end

for i = n:-1:1
    b(i) = (b(i) - LU(i,i+1:n)*b(i+1:n))/LU(i,i);
end

x = b;

end