// Resoluci�n Ax=b, mediante factorizacion LU
N = 4
A=[1 1 0 3; 2 1 -1 1; 3 -1 -1 2; -1 2 3 -1]
B=[4; 1; -3; 4]
X = zeros(N,1);
Y = zeros(N,1);


// Comprobacion
x=inv(A)*B

//B=[8; 7; 14; -7]
// Factorizaci�n LU en la misma matriz A
for i=1:N
   if (A(i,i) == 0.0) then 
     disp("Error en la factorizacion");
     pause;
   end
   for j=i+1:N
      m = A(j,i) / A(i,i);
      A(j, i) = m;
      for k=i+1:N
         A(j, k) = A(j,k) - m * A(i,k);
      end
   end
end
A
B

// Resuelve el sistema LY=B, dados L y B
Y(1) =  B(1) / A(1,1);
for i=2:N
   t = B(i);
   for j=1:i
      t = t - A(i,j) * Y(j);
   end
   Y(i) = t;
end
Y


// Resuelve el sistema UX=Y, dados U  y Y
X(N) = Y(N) / A(N, N);
for i = N-1:-1:1
   t = Y(i);
   for j = N:-1:i
     t = t - A(i,j) * X(j);
   end
   X(i) = t / A(i,i);
end
X