LU factorization without pivoting

The function below (runs in MATLAB and Octave) implements a simple LU factorization. Storing the result in two matrices l and u is mostly for convenience. A regular implementation in a low-level language would simple return the matrix b used in the code.

function [l,u] = simple_lu ( a )
%
%  [l,u] = simple_lu(a) applies the LU factorizaton 
%    to the matrix in a, without pivoting.
%
%  example: a = rand(3,3)
%           [l,u] = simple_lu(a)
%           l*u - a
%
n = size(a,1);
b = a;
for j=1:n-1
  for i=j+1:n
    b(i,j) = b(i,j)/b(j,j);
  end;
  for i=j+1:n
    for k=j+1:n 
      b(i,k) = b(i,k) - b(i,j)*b(j,k);
    end;
  end;
end;
%
% the lower part of the matrix b
% contains l, the upper part contains u.
%
l = zeros(n,n);
u = zeros(n,n);
for i=1:n
  for j=1:n
    if(i>j)
      l(i,j) = b(i,j);
    elseif(i==j)
      l(i,j) = 1;
      u(i,j) = b(i,j);
    else
      u(i,j) = b(i,j);
    end;
  end;
end;