Crout

function [l, u, x]=crout(a, b)

   n=size(a,1)

   a(1,2:n)=a(1,2:n)/a(1,1)

   for m=2:n,

       for k=2:m,

           v=a(m,1:k-1)

           w=a(1:k-1,k)

           a(m,k)=a(m,k)-v*w

       end,

       v=a(m,1:m-1)

       for k=m+1:n,

           w=a(1:m-1,k)

           a(m,k)=(a(m,k)-v*w)/a(m,m)

       end,

   end,

   l=tril(a)

   u=triu(a,1)+eye(n,n)

   y=l\b

   x=u\y

endfunction

function [l, u, x]=crout(a, b)

   n=size(a,1)

   a(1,2:n)=a(1,2:n)/a(1,1)

   for k=2:n,

       for j=2:k

           v=a(k,1:j-1)

           w=a(1:j-1,j)

           a(k,j)=a(k,j)-v*w

       end

       v=a(k,1:k-1)

       for j=k+1:n,

           w=a(1:k-1,j)

           a(k,j)=(a(k,j)-v*w)/a(k,k)

       end

   end

   l=tril(a)

   u=triu(a,1)+eye(n,n)

   y=l\b

   x=u\y

endfunction