function [L,U] = al67(A)
n=length(A);
L=zeros(n,n);
U=zeros(n,n);
U(1,1)=A(1,1);
for i=2:n
 U(i-1,i) = A(i-1, i);
 L(i,i-1) = A(i, i-1)/U(i-1, i-1);
 U(i,i)   = A(i, i)-L(i, i-1)*U(i-1,i);
end
L=L+diag(ones(1,n));

end