clear all;
n=15
A=2*rand(n,n)-1;     % Zufallsmatrix
A=(A+A')/2;      % das macht sie symmetrisch
L=eye(n);        % L muss dimensioniert werden
D=zeros(1,n);    % D muss dimensioniert werden, hier als Zeile
for j=1:n        % Varianten der Formeln von Seite 40 des Buches,
                 % aus A=L*D*L' folgend
    D(1,j)=A(j,j)-(L(j,1:j-1).*D(1,1:j-1))*L(j,1:j-1)';
    for i=j+1:n
        L(i,j)=(A(i,j)-(L(i,1:j-1).*D(1:j-1))*(L(j,1:j-1)'))/D(1,j);   
    end
end
% Testausgabe 
% L
Diagonale=D'
Rest=A-L*diag(D)*L';
restnorm=norm(Rest)