>> % Computation of inv(A) by Gauss-Jordan (see pages 69 and 70)

>> A = [3 2 -1;9 7 -5;-6 6 3]; 
>> A
A =
     3     2    -1
     9     7    -5
    -6     6     3
>> m=[A eye(3)]
m =
     3     2    -1     1     0     0
     9     7    -5     0     1     0
    -6     6     3     0     0     1
>> m(2,:)=m(2,:)-3*m(1,:)
m =
     3     2    -1     1     0     0
     0     1    -2    -3     1     0
    -6     6     3     0     0     1
>> m(3,:)=m(3,:)+2*m(1,:)
m =
     3     2    -1     1     0     0
     0     1    -2    -3     1     0
     0    10     1     2     0     1
>> m(3,:)=m(3,:)-10*m(2,:)
m =
     3     2    -1     1     0     0
     0     1    -2    -3     1     0
     0     0    21    32   -10     1
>> % Quiz: If A=LU what is the matrix on the right? Ans: inv(L)
>> % Let's keep going
>> m(1,:)=m(1,:)-2*m(2,:)
m =
     3     0     3     7    -2     0
     0     1    -2    -3     1     0
     0     0    21    32   -10     1
>> m(1,:)=m(1,:)-(1/7)*m(3,:)
m =
    3.0000         0         0    2.4286   -0.5714   -0.1429
         0    1.0000   -2.0000   -3.0000    1.0000         0
         0         0   21.0000   32.0000  -10.0000    1.0000
>> m(2,:)=m(2,:)+(2/21)*m(3,:)
m =
    3.0000         0         0    2.4286   -0.5714   -0.1429
         0    1.0000         0    0.0476    0.0476    0.0952
         0         0   21.0000   32.0000  -10.0000    1.0000
>> d=diag([3 1 21])
d =
     3     0     0
     0     1     0
     0     0    21
>> inv(d)*m
ans =
    1.0000         0         0    0.8095   -0.1905   -0.0476
         0    1.0000         0    0.0476    0.0476    0.0952
         0         0    1.0000    1.5238   -0.4762    0.0476
>> inv(a)
ans =
    0.8095   -0.1905   -0.0476
    0.0476    0.0476    0.0952
    1.5238   -0.4762    0.0476