Strang

slu.m

function [L,U,P] = slu(A)
[n,n] = size(A); tol = 1.e-6;
P = eye(n);
L = eye(n);
for k = 1 : n
    if abs(A(k,k)) < tol
        r = k + 1;
        while A(r, k) < tol
            r = r + 1;
        end
        A([r k],:) = A([k r],:);
        P([r k],:) = P([k r],:);     
        L([r k],1:k-1) = L([k r],1:k-1);
    end
    L(k,k) = 1;
    for i = k + 1 : n
        L(i,k) = A(i,k) / A(k,k);
        for j = k + 1 : n
            A(i,j) = A(i,j) - L(i,k) * A(k,j);
        end
    end
    for j = k : n
        U(k, j) = A(k, j);
    end
end

slv.m

function x = slv(A,b)
% Solve Ax = b using L and U from slu(A)
[L, U] = slu(A); s = 0;
[n,n] = size(A);
for k = 1 : n
    for j = 1 : k - 1
        s = s + L(k, j) * c(j);
    end
    c(k) = b(k) - s; s = 0;
end
t = 0
for k = n : -1 : 1
    for j = k + 1 : n
        t = t + U(k, j) * x(j);
    end
    x(k) = (c(k) - t) / U(k,k);
end
x = x';