function [L,U]=f_gauss_doolittle(A)
n=length(A);
S=0;
Sb=0;
for k=1:n
  L(k,k)=1;
  for j=k:n
    for s=1:k-1

      S=S+L(k,s)*U(s,j);
    endfor
  U(k,j)=A(k,j)-S;
  endfor
  for i=k+1:n
    for s=1:k-1
      Sb=Sb+(L(i,s)*U(s,k))/U(k,k);
    endfor
  L(i,k)=A(i,k)-Sb;
  endfor
endfor
endfunction


function [x,Ur]=f_gauss(A,b)
Ur=[A,b];
n=length(b);
  for k=1:1:n-1;
    if A(k,k)==0
      error('Blad');
    endif
    for j=k+1:1:n+1
      for i=k+1:1:n
      Ur(i,j)=Ur(i,j)-(Ur(i,k))/(Ur(k,k))*Ur(k,j);
      endfor
    endfor
  endfor
x=0;
  for i=n:-1:1;
    S=0;
    for j=i+1:1:n;
      S=S+Ur(i,j)*x(j);
  endfor
  x(i)=(Ur(i,n+1)-S)/(Ur(i,i));
  endfor
endfunction


function [x,Ur]=f_gauss(A,b)
Ur=[A,b];
n=length(b);
  for k=1:1:n-1;
    if A(k,k)==0
      error('Blad');
    endif
    for j=k+1:1:n+1
      for i=k+1:1:n
      Ur(i,j)=Ur(i,j)-(Ur(i,k))/(Ur(k,k))*Ur(k,j);
      endfor
    endfor
  endfor
x=0;
  for i=n:-1:1;
    S=0;
    for j=i+1:1:n;
      S=S+Ur(i,j)*x(j);
  endfor
  x(i)=(Ur(i,n+1)-S)/(Ur(i,i));
  endfor
endfunction


A=[3,1,-1;-1,5,-1;2,4,8]
b=[6;10;2]
A\b
[x,Ur]=f_gauss(A,b)
[L,U]=lu(A)
[L,U]=f_gauss_doolittle(A)