Método de LU 6x6

les dejo un código hecho en MATLAB del método de resolución para ecuaciones lineales en un sistema de 6 incógnitas acomodado en una matriz principal de 6x6 y sus términos independientes en un vector b[6], bien lo pude hacer para nxn pero solo me es útil por este momento este tamaño

%METODO DE LU

%EMMANUEL (YETLA)

%ANALISIS NUMERICO

clear,clc

fprintf(' Metodo de LU');

fprintf('\nDe la ecuacion numero 1 ingresa los coeficientes');

for n=1:6;

fprintf('\n x%d= ',n);

r1(n)=input('');

end;

b(1)=input('\n b= ');

fprintf('De la ecuacion numero 2 ingresa los coeficientes');

for n=1:6;

fprintf('\n x%d= ',n);

r2(n)=input('');

end;

b(2)=input('\n b= ');

fprintf('De la ecuacion numero 3 ingresa los coeficientes');

for n=1:6;

fprintf('\n x%d= ',n);

r3(n)=input('');

end;

b(3)=input('\n b= ');

fprintf('De la ecuacion numero 4 ingresa los coeficientes');

for n=1:6;

fprintf('\n x%d= ',n);

r4(n)=input('');

end;

b(4)=input('\n b= ');

fprintf('De la ecuacion numero 5 ingresa los coeficientes');

for n=1:6;

fprintf('\n x%d= ',n);

r5(n)=input('');

end;

b(5)=input('\n b= ');

fprintf('De la ecuacion numero 6 ingresa los coeficientes');

for n=1:6;

fprintf('\n x%d= ',n);

r6(n)=input('');

end;

b(6)=input('\n b= ');

fprintf('\nMatriz A');

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d \t = %d',r1(1),r1(2),r1(3),r1(4),r1(5),r1(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d \t = %d',r2(1),r2(2),r2(3),r2(4),r2(5),r2(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d \t = %d',r3(1),r3(2),r3(3),r3(4),r3(5),r3(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d \t = %d',r4(1),r4(2),r4(3),r4(4),r4(5),r4(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d \t = %d',r5(1),r5(2),r5(3),r5(4),r5(5),r5(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d \t = %d',r6(1),r6(2),r6(3),r6(4),r6(5),r6(6));

fprintf('\nMatriz X');

fprintf('\n F1');

fprintf('\n F2');

fprintf('\n F3');

fprintf('\n H2');

fprintf('\n V2');

fprintf('\n V3');

fprintf('\nMatriz b');

fprintf('\n %d ',b(1));

fprintf('\n %d ',b(2));

fprintf('\n %d ',b(3));

fprintf('\n %d ',b(4));

fprintf('\n %d ',b(5));

fprintf('\n %d ',b(6));

a1=r1;

%columna 1

f1=r2(1)/r1(1);

a2=r2-r1*f1;

f2=r3(1)/r1(1);

a3=r3-r1*f2;

f3=r4(1)/r1(1);

a4=r4-r1*f3;

f4=r5(1)/r1(1);

a5=r5-r1*f4;

f5=r6(1)/r1(1);

a6=r6-r1*f5;

%columna 2

f6=a3(2)/a2(2);

a3=a3-a2*f6;

f7=a4(2)/a2(2);

a4=a4-a2*f7;

f8=a5(2)/a2(2);

a5=a5-a2*f8;

f9=a6(2)/a2(2);

a6=a6-a2*f9;

%columna 3

f10=a4(3)/a3(3);

a4=a4-a3*f10;

f11=a5(3)/a3(3);

a5=a5-a3*f11;

f12=a6(3)/a3(3);

a6=a6-a3*f12;

%columna 4

f13=a5(4)/a4(4);

a5=a5-a4*f13;

f14=a6(4)/a4(4);

a6=a6-a4*f14;

%columna 5

f15=a6(5)/a5(5);

a6=a6-a5*f15;

a2(1)=0;

a3(1)=0;a3(2)=0;

a4(1)=0;a4(2)=0;a4(3)=0;

a5(1)=0;a5(2)=0;a5(3)=0;a5(4)=0;

a6(1)=0;a6(2)=0;a6(3)=0;a6(4)=0;a6(5)=0;

fprintf('\nMatriz U');

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',a1(1),a1(2),a1(3),a1(4),a1(5),a1(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',a2(1),a2(2),a2(3),a2(4),a2(5),a2(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',a3(1),a3(2),a3(3),a3(4),a3(5),a3(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',a4(1),a4(2),a4(3),a4(4),a4(5),a4(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',a5(1),a5(2),a5(3),a5(4),a5(5),a5(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',a6(1),a6(2),a6(3),a6(4),a6(5),a6(6));

b1=[1 0 0 0 0 0];

b2=[f1 1 0 0 0 0];

b3=[f2 f6 1 0 0 0];

b4=[f3 f7 f10 1 0 0];

b5=[f4 f8 f11 f13 1 0];

b6=[f5 f9 f12 f14 f15 1];

fprintf('\nMatriz L');

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',b1(1),b1(2),b1(3),b1(4),b1(5),b1(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',b2(1),b2(2),b2(3),b2(4),b2(5),b2(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',b3(1),b3(2),b3(3),b3(4),b3(5),b3(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',b4(1),b4(2),b4(3),b4(4),b4(5),b4(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',b5(1),b5(2),b5(3),b5(4),b5(5),b5(6));

fprintf('\n %d \t %d \t %d \t %d \t %d \t %d',b6(1),b6(2),b6(3),b6(4),b6(5),b6(6));

c1=b(1);

c2=b(2)-c1*f1;

c3=b(3)-c1*f2-c2*f6;

c4=b(4)-c1*f3-c2*f7-c3*f10;

c5=b(5)-c1*f4-c2*f8-c3*f11-c4*f13;

c6=b(6)-c1*f5-c2*f9-c3*f12-c4*f14-c5*f15;

fprintf('\n c1 = %d ',c1);

fprintf('\n c2 = %d ',c2);

fprintf('\n c3 = %d ',c3);

fprintf('\n c4 = %d ',c4);

fprintf('\n c5 = %d ',c5);

fprintf('\n c6 = %d ',c6);

x6=c6/a6(6)

x5=(c5-x6*a5(6))/a5(5)

x4=(c4-x6*a4(6)-x5*a4(5))/a4(4)

x3=(c3-x6*a3(6)-x5*a3(5)-x4*a1(4))/a3(3)

x2=(c2-x6*a2(6)-x5*a2(5)-x4*a2(4)-x3*a2(3))/a2(2)

x1=(c1-x6*a1(6)-x5*a1(5)-x4*a1(4)-x3*a1(3)-x2*a1(2))/a1(1)