function x=lusol(A,b)
[L,U,P]=lu(A);
d=P*b;
n = length(b);
y = zeros(n,1);
for j=1:n-1
   y(j) = d(j)/L(j,j);
   d(j+1:n) = d(j+1:n) - L(j+1:n,j)*y(j);
end
y(n) = d(n)/L(n,n);

x = zeros(n,1);
for j=n:-1:2
   x(j) = y(j)/U(j,j);
   y(1:j-1) = y(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = y(1)/U(1,1);