>> a=[1 3 0 3]

a =

     1     3     0     3


>> A1=[-1 -1 -1 1; 0 -4 2 -8; 2 0 3 -1] 

A1 =

    -1    -1    -1     1
     0    -4     2    -8
     2     0     3    -1

> A =[ a;A1]

A =

     1     3     0     3
    -1    -1    -1     1
     0    -4     2    -8
     2     0     3    -1

>> I4 = eye(4)

I4 =

     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

>> e1=I4(1,1:4)

e1 =

     1     0     0     0

>> e2=I4(2,1:4)

e2 =

     0     1     0     0

>> e3=I4(3,1:4)

e3 =

     0     0     1     0

>> e4=I4(4,1:4)

e4 =

     0     0     0     1


>> R41_2=[e1;e2;e3;e4-2*e1]

R41_2 =

     1     0     0     0
     0     1     0     0
     0     0     1     0
    -2     0     0     1

>> A=R41_2*A

A =

     1     3     0     3
    -1    -1    -1     1
     0    -4     2    -8
     0    -6     3    -7

R211 =

     1     0     0     0
     1     1     0     0
     0     0     1     0
     0     0     0     1

>> A=R211*A

A =

     1     3     0     3
     0     2    -1     4
     0    -4     2    -8
     0    -6     3    -7


R322 =

     1     0     0     0
     0     1     0     0
     0     2     1     0
     0     0     0     1

>> A=R322*A

A =

     1     3     0     3
     0     2    -1     4
     0     0     0     0
     0    -6     3    -7

>> Sw34=[e1;e2;e4;e3]

Sw34 =

     1     0     0     0
     0     1     0     0
     0     0     0     1
     0     0     1     0

>> A=Sw34*A

A =

     1     3     0     3
     0     2    -1     4
     0    -6     3    -7
     0     0     0     0

>> R323=[e1;e2;e3+3*e2;e4]

R323 =

     1     0     0     0
     0     1     0     0
     0     3     1     0
     0     0     0     1

>> A=R323*A

A =

     1     3     0     3
     0     2    -1     4
     0     0     0     5
     0     0     0     0

Forward phase is complete.  A is in Row Echelon Form. Sometimes this is ref(A).  Note: Scale(R1,c) was unnecessary for this phase.

Now for the Backward Phase.

>> Sc2d2=[e1;e2*0.5;e3;e4]

Sc2d2 =

    1.0000         0         0         0
         0    0.5000         0         0
         0         0    1.0000         0
         0         0         0    1.0000

>> A=Sc2d2*A

A =

    1.0000    3.0000         0    3.0000
         0    1.0000   -0.5000    2.0000
         0         0         0    5.0000
         0         0         0         0

>> R12_3=[e1-3*e2;e2;e3;e4]

R12_3 =

     1    -3     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

>> A=R12_3*A

A =

    1.0000         0    1.5000   -3.0000
         0    1.0000   -0.5000    2.0000
         0         0         0    5.0000
         0         0         0         0

Sc3d5 =

    1.0000         0         0         0
         0    1.0000         0         0
         0         0    0.2000         0
         0         0         0    1.0000

>> A=Sc3d5*A

A =

    1.0000         0    1.5000   -3.0000
         0    1.0000   -0.5000    2.0000
         0         0         0    1.0000
         0         0         0         0

NOTE: Here I named the elementary matrix incorrectly.

>> R133=[e1+3*e3;e2;e3;e4]

R133 =

     1     0     3     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

>> A=R133*A

A =

    1.0000         0    1.5000         0
         0    1.0000   -0.5000    2.0000
         0         0         0    1.0000
         0         0         0         0

>> R23_2=[e1;e2-2*e3;e3;e4]

R23_2 =

     1     0     0     0
     0     1    -2     0
     0     0     1     0
     0     0     0     1

>> A=R23_2*A

A =

    1.0000         0    1.5000         0
         0    1.0000   -0.5000         0
         0         0         0    1.0000
         0         0         0         0


A is in reduced row echelon form the backward phase is complete.