define a random 5x5 matrix:
>> A = rand(5,5);
use _my_ lu decomposition program to find the lower and upper
triangular matrices L and U.
>> [L,U] = my_lu(A);
Here they are:
>> L
L =
    1.0000         0         0         0         0
    0.9801    1.0000         0         0         0
    2.9778    0.5461    1.0000         0         0
    1.3424    1.2453    0.9402    1.0000         0
    0.9805    0.6163    0.0073    0.3323    1.0000
>> U
U =
    0.2028    0.0153    0.4186    0.8381    0.5028
         0    0.7318    0.4359   -0.8018    0.2167
         0         0   -0.9595   -1.3767   -1.1867
         0         0         0    1.5472    0.4756
         0         0         0         0   -0.5863

Test that we did it right: that L*U = A:
>> max(max(abs(L*U-A)))
ans =
   3.3307e-16

Now we use matlab's lu decomposition and check that it's right...
>> [L,U] = lu(A);
>> max(max(abs(L*U-A)))
ans =
   1.1102e-16

We look at the matrices:
>> L, U
L =
    0.3358   -0.1835    0.3365    0.9247    1.0000
    0.3291    0.8210    1.0000         0         0
    1.0000         0         0         0         0
    0.4508    1.0000         0         0         0
    0.3293    0.4369    0.7330    1.0000         0
U =
    0.6038    0.4451    0.5252    0.6813    0.4289
         0    0.7312   -0.0341    0.0724    0.1113
         0         0    0.7014   -0.2640    0.4770
         0         0         0    0.7694   -0.3498
         0         0         0         0    0.5422

Whoah!  That L isn't lower triangular.  Look at what matlab does:
>> help lu
 LU     LU factorization.
    [L,U] = lu(X) stores an upper triangular matrix in U and a
    "psychologically lower triangular matrix" (i.e. a product
    of lower triangular and permutation matrices) in L, so
    that X = L*U.  X must be square.
 
    [L,U,P] = lu(X) returns lower triangular matrix L, upper
    triangular matrix U, and permutation matrix P so that
    P*X = L*U.
 
    LU(X), with one output argument, returns the output from
    LINPACK'S ZGEFA routine.
 
    LU(X,THRESH) controls pivoting in sparse matrices, where
    THRESH is a pivot threshold in [0,1].  Pivoting occurs
    when the diagonal entry in a column has magnitude less than
    THRESH times the magnitude of any sub-diagonal entry in that
    column.  THRESH = 0 forces diagonal pivoting.  THRESH = 1 is
    the default.
 
    See also LUINC, QR, RREF.

It refers to a "psychologically lower triangular matrix".  If we
ask for three outputs then we get the lower triangular and upper
triangular matrices as well as the permutation matrix P.  But
now we don't have L*U = A, we have L*U = P*A:
>> [L,U,P] = lu(A)   
L =
    1.0000         0         0         0         0
    0.4508    1.0000         0         0         0
    0.3291    0.8210    1.0000         0         0
    0.3293    0.4369    0.7330    1.0000         0
    0.3358   -0.1835    0.3365    0.9247    1.0000
U =
    0.6038    0.4451    0.5252    0.6813    0.4289
         0    0.7312   -0.0341    0.0724    0.1113
         0         0    0.7014   -0.2640    0.4770
         0         0         0    0.7694   -0.3498
         0         0         0         0    0.5422
P =
     0     0     1     0     0
     0     0     0     1     0
     0     1     0     0     0
     0     0     0     0     1
     1     0     0     0     0

>> max(max(abs(L*U - P*A)))
ans =
   1.1102e-16

We can check that their LU decomposition is better than mine...
>> n = 100; A = rand(n,n);
>> [L,U] = my_lu(A); max(max(abs(L*U-A)))
ans =
   1.9784e-13
>> [L,U] = lu(A); max(max(abs(L*U-A)))
ans =
   1.8319e-15


>> diary off