Crout's algorithm works just fine

Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m

@@ -1,9 +1,9 @@
 function A = Alg1_outer_product_gaussian_elimination(A)
 % ADDME Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1)
+% Matrix A has to be a square singular matrix.
 
-
+[m, n] = size(A);
-[n, m] = size(A);
+if m ~= n
-if n ~= m
     error('Matrix is not squared!')
 end
     
@@ -11,8 +11,8 @@ end
 %     error('Matrix is not nonsingular!')
 %  end
  
- for k = 1 : n-1
+ for k = 1 : m-1
-     rows = k + 1 : n;
+     rows = k + 1 : m;
      A(rows, k) = A(rows, k)/A(k, k);
      A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
  end

Direct Methods for Solving Linear Systems/Alg7.m

@@ -5,9 +5,30 @@ function [L, U] = Alg7(A)
 
 [m, n] = size(A);
 
-A = Alg1_outer_product_gaussian_elimination(A);
+% This solution is mathematically correct, however computationaly
+% inefficient.
+%
+% A = Alg1_outer_product_gaussian_elimination(A);
+% 
+% U = triu(A);
+% L = tril(A, -1) + eye(m);
 
-U = triu(A);
+% Instead, we should use the Crout's algorithm
-L = tril(A, -1) + eye(m);
+% It returns an Unit Upper Triangular matrix and a Lower Traingular matrix
+% (in oppose to the Gaussian Elimination).
+A
+
+L = zeros(m, n)
+U = eye(m, n)
+
+for k = 1:n
+    for i = k:n
+        L(i, k) = A(i, k) - dot(L(i,1:k-1), U(1:k-1, k));
+    end
+    
+    for j = k : n
+        U(k, j) = (A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j))) / L(k, k);
+    end
+end
 
 end