function [x,L,U,error]=gauss(A,b)
N=length(A);
x=zeros(N,1);
L=zeros(N);
U=zeros(N);
error=0;
A=[A,b];
for k=1:N-1
    if A(k,k)==0
        error=1;
        for i=k+1:N
            if A(i,k)~=0, aux=A(i,:); A(i,:)=A(k,:); A(k,:)=aux; error=0; break, end
        end
        if error==1, return, end
    end
    i=[k+1:N];
    j=[k+1:N+1];
    A(i,k)=A(i,k)/A(k,k);
    A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
if A(N,N)==0, error=1; return, end
c=A(:,N+1);
A=A(:,1:N);
U=triu(A);
L=eye(N)+tril(A,-1);
x(N)=c(N)/U(N,N);
for i=N-1:-1:1
    j=[i+1:N];
    x(i)=(c(i)-U(i,j)*x(j))/U(i,i);
end