clear; clc;

told=1.0;
for n = [ 50 100 200 400 800 1600 3200  6400  12800 ];
 
  a=rand(n,n); 
  b=rand(n,1);

  tic;  %  t0=cputime;
    xx=mylin(a,b,n);
  %  xx=a\b;
  %  t=cputime-t0;
  t=toc;
  ratio=t/told;
  fprintf(' n=%6d   time= %12.4f  ratio=%12.4f\n',n,t,ratio);
  told=t;
  if t>10; break; end;
end

function x=mylin(a,b,n)
for ic=1:(n-1)   %..... loop over all columns except last
    for ir=(ic+1):n        %...... loop over rows below the diagonal
        mult=a(ir,ic)/a(ic,ic);
        for icc=ic:n
            a(ir,icc) = a(ir,icc) - mult*a(ic,icc);  %  subtract multiple of above row from this row
        end
        b(ir) = b(ir) - mult*b(ic);  %  do the same with the right hand side
    end
end
%  a should be all zeros below the diagonal now

 %.... back substitution
x(n)=b(n)/a(n,n);
for ir=(n-1):-1:1  %   backwards over rows
    for ic=(ir+1):n
        b(ir) = b(ir) - a(ir,ic)*x(ic); % subtract the known part
    end   
    x(ir)=b(ir)/a(ir,ir);  %   and divide by the diagonal
end

end