function [x,y]=Crout(A,b)

n=length(A(1,:));
if length(A(:,1))~=n
    error('Macierz nie jest kwadratowa!');
    return;
end
if length(b)~=n
    error('Podany wektor ma nieprawidłowy rozmiar!');
    return;
end

U=eye(n);

L=zeros(n);
y=diag(L);
x=diag(L);


%obliczanie L i U

for k=1:n
    for i=k:n

        sumL=0.0;
        for p=1:(k-1)
            sumL=sumL+L(i,p).*U(p,k);
        end

        L(i,k)=A(i,k)-sumL;
    end

    for j=(k+1):n

        sumU=0.0;
        for p=1:(k-1)
            sumU=sumU+L(k,p).*U(p,j);
        end

        U(k,j)=(A(k,j)-sumU)./L(k,k);
    end
end

%obliczanie rozwiązania

y(1)=b(1);

for i=2:n

    sumY=0.0;
    for k=1:i-1
        sumY=sumY+L(i,k).*y(k);
    end

    y(i)=b(i)-sumY;
end

x(n)=y(n)/U(n,n);

for i=n-1:-1:1
    sumX=0.00;
    for k=(i+1):n
        sumX=sumX+U(i,k).*x(k);
    end
    x(i)=(y(i)-sumX)./U(i,i);
end

end