% random equations
ndim = 5
A = rand(ndim);
b = rand(ndim,1);

% solution using matlab's '\' operator
z = A\b;

% solution using 'homemade' lu decomposition
[L, U, p] = ludecomp(A);
y = forward(L, b(p()));
x = backsubs(U, y);

% How good is the agreement?
max(abs(z-x))
