function [ x ] = Asolve( A,b )
%LUsolve: solve Ax=b by LU decomposition
  [L,U] = LUdecomp(A);
  y = Lsolve(L,b);
  x = Usolve(U,y);
end

