function [B, D, C] = crout(A)
n=length(A);
B=zeros(1,n);    % sub-diagonal of L. B(i) will be L(i,i-1), i = 2, ..., n. 
                 % B(1) is dummy
D=zeros(1,n);    % diaginal of U. D(i) will be U(i,i), i = 1, ..., n.
C=zeros(1,n-1);  % super-diagonal of U. C(i) will be U(i,i+1), i = 1, ..., n-1. 
                 % See also problem 27 on page 413.
D(1) = A(1,1);   % U(1,1) = A(1,1)
for i=2:n
 C(i-1) = A(i-1,i);             % U(i-1,i) = A(i-1,i);
 B(i)   = A(i,i-1)/D(i-1);      % L(i,i-1) = A(i,i-1)/U(i-1,i-1);
 D(i)   = A(i,i)-B(i)*C(i-1);   % U(i,i)   = A(i,i)-L(i,i-1)*U(i-1,i);
end

end