Doolitle
function [x]=doolittle(a,b)
  n=size(a,1)
  l=eye(n,n)
  u=zeros(n,n)

  for i=1:n,
    u(1,i)=a(1,i);
  end
  for i=2:n,
    l(i,1)=a(i,1)/u(1,1);
  end
  for j=2:n-1,
    for k=j:n,
      y=0;
      for s=1:j-1,
        y=y+l(j,s)*u(s,k);
      end
      u(j,k)=a(j,k)-y;
      for q=j:n-1,
        z=0;
        for p=q+1:n,
          z=z+l(p,s)*u(s,q);
          l(p,q)=1/u(q,q)*(a(p,q)-z);
        end
      end
    end
  end
  w=0;
  for s=1:n-1,
    w=w+l(n,s)*u(s,n);
  end
  u(n,n)=a(n,n)-w;
  y=l\b;
  x=u\y;
endfunction

Crout
function [x]=crout(a,b)
  n=size(a,1)
  l=zeros(n,n)
  u=eye(n,n)

  for i=1:n,
    l(i,1)=a(i,1);
  end
  for i=2:n,
    u(1,i)=a(1,i)/l(1,1);
  end
  for k=2:n-1,
    for j=k:n,
      y=0;
      for s=1:k-1,
        y=y+l(j,s)*u(s,k);
      end
      l(j,k)=a(j,k)-y;
      for p=k:n-1,
        z=0;
        for q=p+1:n,
          z=z+l(p,s)*u(s,q);
          u(p,q)=1/l(p,p)*(a(p,q)-z);
        end
      end
    end
  end
  w=0;
  for s=1:n-1,
    w=w+l(n,s)*u(s,n);
  end
  l(n,n)=a(n,n)-w;
  y=l\b;
  x=u\y;
endfunction

Hasil Eksekusi
-->A=[3 -1 1;10 3 -2;1 1 -1]
 A  =
    3.   - 1.    1.
    10.    3.  - 2.
    1.     1.  - 1.

-->A=[3 -1 1;10 3 -2;1 1 -1];

-->B=[12;-23;6];

-->getd

-->x=doolittle(A,B)
 x  =
    4.5
  - 71.
  - 72.5

-->x=crout(A,B)
 x  =
    4.5
  - 71.
  - 72.5

-->x=A\B
 x  =
    4.5
  - 71.
  - 72.5

-->A=[2 6 -1;1 -8 3;1 6 -10];

-->B=[0;0;0];

-->getd
Warning :redefining function: crout
Warning :redefining function: doolittle

-->x=doolittle(A,B)
 x  =
    0.
    0.
    0.

-->x=crout(A,B)
 x  =
    0.
    0.
    0.

-->x=A\B
 x  =
    0.
    0.
    0.

-->A=[1 1 -3;2 2 -6;1 -9 2];

-->B=[9;18;3];

-->getd
Warning :redefining function: crout
Warning :redefining function: doolittle

-->x=doolittle(A,B)
 !--error 27
division by zero...
at line      23 of function doolittle called by :
x=doolittle(A,B)

-->x=crout(A,B)
 !--error 27
division by zero...
at line      23 of function crout called by :
x=crout(A,B)

-->x=A\B
 warning
 matrix is close to singular or badly scaled. rcond =    0.0000D+00
computing least squares solution. (see lsq)

 x  =
    0.
  - 1.08
  - 3.36

-->diary(0)