function [L,U] = my_lu( eps )

%   [ eps -1;
%      1   1 ]
 
l21 = 1/single(eps);

u11 = eps;
u12 = -1;
u22 = 1 + 1/single(eps);

L = [1   0; 
     l21 1 ];

U = [u11 u12;
     0   u22];

end



%  try
%  eps = 10^-7; A = [eps -1; 1 1]; [L,U] = my_lu(eps); norm(L*U-A)
%  eps = 10^-8; A = [eps -1; 1 1]; [L,U] = my_lu(eps); norm(L*U-A)
