%risoluzione LU con Doolittle e Crout

clear all
clc

A=[4 1 2; 1 3 1; 0 1 2];
b=sum(A,2);

[L,U] = DOOLITTLE(A);
[L,U]=CROUT(A);


function [L,U] = DOOLITTLE(A)
n=length(A);
L=eye(n);
U(1,:)=A(1,:);

for k=1:n
    U(k,k:n)=A(k,k:n) -L(k,1:k-1)*U(1:k-1,k:n);
    L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
end

end

function [L,U] = CROUT(A)
n=length(A);
U=eye(n);
L(:,1)=A(:,1);

for k=1:n
    L(k:n,k)=A(k:n,k) -L(k:n,1:k-1)*U(1:k-1,k);
    U(k,k+1:n)=(A(k,k+1:n)-L(k,1:k-1)*U(1:k-1,k+1:n))/L(k,k);
end
end