A=[67 -13 -28 -26                  % example of Gauss 
   -13 69 -50 -6
   -28 -50 156 -78
   -26 -6 -78 110];
f=-[6                              % original right hand side
  -7558
  -14604
   22156];
id=[4 2 1 3 1 2 1];                % sequence of variables used
D=diag(A);
fn=f;
u=zeros(4,1);                      % to sum corrections
F=[];
for i=1:length(id)                 % compute correction steps
  d=zeros(4,1);
  d(id(i))=round(fn(id(i))/D(id(i)))
  u(id(i))=u(id(i))+d(id(i));      % add correction
  fn=fn-A*d;                       % update the right-hand side
  F=[F -fn]                        % collect updated right-hand sides
  pause
end
u                                  % approximate solution of Gauss
uexact=A\f                         % compare with Matlab solution
[A*u-f A*uexact-f]                 % Gauss and Matlab solution residual
fm=f-u(4)*A(:,4);                  % compute Matlab solution with same 4th 
uc=A(1:3,:)'\fm;                   % component as Gauss' solution
uc(4)=u(4);
[u uc]                             % compare again
