// define matrix A for system
// x+y=5
// 2x+y=4
A=[1 1 5; 2 1 4]
// row 1
A(1,:)
// row 2
A(2,:)
// col 1
A(:,1)
// entry 1,3 = entry in row 1, col 3
A(1,3)
// row reduce A
A(2,:)=A(2,:)-2*A(1,:)
A(2,:)=-1*A(2,:)
// this corresponds to the equivalent system
// x+y=5
//   y=6
// solve using back substitution
y=6
x=5-y



//
// define matrix B for system
// x+y+z=0
// x+2y-z=5
// -2x+6y+4z=2
B=[1 1 1 0; 1 2 -1 5; -2 6 4 2]
// row reduce B
B(2,:)=B(2,:)-B(1,:)
B(3,:)=B(3,:)+2*B(1,:)
B(3,:)=B(3,:)-8*B(2,:)
B(3,:)=1/22*B(3,:)
// solve using back sub
z=B(3,4)
y=B(2,4)+2*z
x=B(1,4)-y-z
// check that it solves original system
x+y+z
x+2*y-z
-2*x+6*y+4*z



//
// system with no solutions
// 2x+3y=10
// -2x-3y=5
C=[2 3 10; -2 -3 5]
// row reduce C
C(2,:)=C(1,:)+C(2,:)
// equivalent system is
// 2x+3y=10
// 0x+0y=15
// the last equation says 0=15, so there are no solutions



//
// system with infinitely many solutions
// x+y-z=10
// -2x+3y+3z=11
// -3x+7y+3z=32
// define matrix D for system
D=[1 1 -1 10; -2 3 3 11; -3 7 5 32]
// row reduce
D(2,:)=D(2,:)+2*D(1,:)
D(3,:)=D(3,:)+3*D(1,:)
D(3,:)=D(3,:)-2*D(2,:)
D(2,:)=D(2,:)/5
// new system is
// x+y-z=10
//   y+.2z=6.2
//       0=0
// solve using back sub
// z is free
// introduce a parameter t, and let z=t
// x, y, z then become functions of t
function [x, y ,z]=solution(t)
	 z=t
	 y=6.2-.2*z
	 x=10-y+z
endfunction
// different values of t will give different solutions
[x y z]=solution(1)
[x y z]=solution(-2)
// as a column vector (; at end of line suppresses output of a command)
[x y z]=solution(4.1);
[x;y;z]


//
// more complicated example
// 6 unknowns x1 x2 x3 x4 x5 x6, 4 equations
// represented by a 6 by 4 matrix
E= [ 1 1 1 0 -1 2 3;
     5 5 -10 2.5 1 -4 3.1;
     2 3 4 5 6 7 8;
     0 -2.1 3.7 6 2 4 5 ]
// zero out first column in rows 2,3,4
E(2,:)=E(2,:)-5*E(1,:)
E(3,:)=E(3,:)-2*E(1,:)
// switch rows 2 and 3
E([2 3],:)=E([3 2],:)
// zero out second column in rows 3, 4
E(4,:)=E(4,:)+2.1*E(2,:)
// zero out third column in row 4
E(4,:)=E(4,:)*15+7.9*E(3,:)
// normalize leading entries to 1
E(3,:)=E(3,:)/E(3,3)
E(4,:)=E(4,:)/E(4,4)
// back sub to find solution
// x5 x6 are free variables, call them s and t
function [x1, x2, x3, x4, x5, x6]=solution(s,t)
	 x5=s
	 x6=t
	 x4=E(4,7)-t*E(4,6)-s*E(4,5)
	 x3=E(3,7)-t*E(3,6)-s*E(3,5)-x4*E(3,4)
	 x2=E(2,7)-t*E(2,6)-s*E(2,5)-x4*E(2,4)-x3*E(2,3)
	 x1=E(1,7)-t*E(1,6)-s*E(1,5)-x4*E(1,4)-x3*E(1,3)-x2*E(1,2)
endfunction
// some sample solutions
[x1,x2,x3,x4,x5,x6]=solution(1,0);
[x1;x2;x3;x4;x5;x6]
[x1,x2,x3,x4,x5,x6]=solution(-5,11);
[x1;x2;x3;x4;x5;x6]
[x1,x2,x3,x4,x5,x6]=solution(3.14,-2.5);
[x1;x2;x3;x4;x5;x6]