n=input('n=')

 a= rand(n,n); xe=ones(n,1); b=a*xe;


 [L, U, x] = my_LU_solver(n, a, b)

% Check LU 

 norm(a-L*U)

% check solution

 norm(x-xe)