diff --git a/Direct Methods for Solving Linear Systems/Alg11.m b/Direct Methods for Solving Linear Systems/Alg11.m
index 3fd774b..b33c7c4 100644
--- a/Direct Methods for Solving Linear Systems/Alg11.m	
+++ b/Direct Methods for Solving Linear Systems/Alg11.m	
@@ -1,5 +1,5 @@
 function [Q R] = Alg11(A)
-% Algorithm 11: QR factorization via Housholder algorithm.
+% Algorithm 11: QR factorization via Householder algorithm.
 
 [m, n] = size(A);
 
@@ -17,5 +17,8 @@ end
 
 R = triu(A)
 
+H = tril(A, -1)
+
+
 end
 
diff --git a/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m b/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m
index 31648c1..babce75 100644
--- a/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m	
+++ b/Direct Methods for Solving Linear Systems/Alg1_outer_product_gaussian_elimination.m	
@@ -1,21 +1,25 @@
 function A = Alg1_outer_product_gaussian_elimination(A)
-% Algorithm 1: Outer Product Gaussian Elimination (Golub, Loan, 3.2.1)
+% Algorithm 1: Outer Product Gaussian Elimination
 % Performs a gaussian eliminaion on a square matrix A.
 
-    [m, n] = size(A);
+[m, n] = size(A);
 
-    if m ~= n
-        error('Matrix is not squared!')
-    end
-    
-%  if det(A) == 0
-%     error('Matrix is not nonsingular!')
-%  end
- 
-    for k = 1 : m-1
-     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);
+if m ~= n
+    error('Matrix is not square!')
+end
+
+for k = 1:n-1
+    if det(A(1:k, 1:k)) < eps
+        error('Matrix is not nonsingular!')
     end
+end
+
+% The following algorithm is based on the Algrotihm 3.2.1 from [2].
+
+for k = 1 : m-1
+    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
  
 end
\ No newline at end of file
diff --git a/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m b/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m
index 7c01c8e..43f9cfe 100644
--- a/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m	
+++ b/Direct Methods for Solving Linear Systems/Alg2_gaussian_elimination_with_complete_pivoting.m	
@@ -1,43 +1,47 @@
 function [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A)
 % Algorithm 2: Gaussian Elimination with Complete Pivoting.
-% [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A) 
+% [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A)
+% computes the complete pivoting factorization PAQ = LU.
 
-[n, m] = size(A);
+[m, n] = size(A);
 
-if n ~= m
-    error('Matrix is not squared!')
+if m ~= n
+    error('Matrix is not square!')
 end
-    
-% if det(A) == 0
-%     error('Matrix is not nonsingular!')
-% end
 
-p = 1:n;
-q = 1:n;
+% p and q are permutation vectors – respectively rows and columns
+p = 1:m;
+q = 1:m;
 
-for k = 1 : n-1
-    i = k:n;
-    j = k:n;
+% The following algorithm is based on the Algrotihm 3.4.2 from [2].
+
+for k = 1 : m-1
+    i = k:m;
+    j = k:m;
+    % Find the maximum entry to be the next pivot
     [max_val, rows_of_max_in_col] = max(abs(A(i, j)));
     [~, max_col] = max(max_val);
     max_row = rows_of_max_in_col(max_col);
-    % Assign value of mu and lambda in respect to the main A matrix
+    % Assign value of mu and lambda in respect to the main matrix A
     [mi, lm] = deal(max_row+k-1, max_col+k-1);
-    A([k mi], 1:n) = deal(A([mi k], 1:n));
-    A(1:n, [k lm]) = deal(A(1:n, [lm k]));
+    % Interchange the rows and columns of matrix A...
+    A([k mi], 1:m) = deal(A([mi k], 1:m));
+    A(1:m, [k lm]) = deal(A(1:m, [lm k]));
+    % ...and respective permutation vectors entries.
     p([k, mi]) = p([mi, k]);
     q([k, lm]) = q([lm, k]);
-    
     % Perform Gaussian elimination with the greatest pivot
     if A(k, k) ~= 0
-        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
 end
 
-I = eye(n);
+I = eye(m);
 U = triu(A);
 L = tril(A, -1) + I;
 P = I(p, :);
-Q = I(:, q);
\ No newline at end of file
+Q = I(:, q);
+
+end
\ No newline at end of file
diff --git a/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m b/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m
index a0618c0..b793584 100644
--- a/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m	
+++ b/Direct Methods for Solving Linear Systems/Alg3_forward_substitution.m	
@@ -1,18 +1,29 @@
 function b = Alg3_forward_substitution(L, b)
-% Algorithm 3: Forward Substitution (Golub, Loan, Alg. 3.1.1)
+% Algorithm 3: Forward Substitution
+% b = Alg3_forward_substitution(L, b) overwrites b with the solution to
+% Lx = b.
 
-    [m, n] = size(L);
-    if m ~= n
-        error('Matrix is not squar!')
-    end
+[m, n] = size(L);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+if length(b) ~= m
+    error('Vector b has wrong length!')
+end
+
+if det(L) < eps
+    error("Matrix is not nonsingular!")
+end
+
+% The following algorithm is based on the Algrotihm 3.1.1 from [2].
+
+% b(m, :) so that matrices are also accepted
+b(1, :) = b(1, :)/L(1,1);
+for i = 2:m
+    b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i);
+end
 
-    if length(b) ~= m
-        error('Vector b has wrong length!')
-    end
-
-    b(1, :) = b(1, :)/L(1,1);
-    for i = 2:m
-        b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i);
-    end
 end
 
diff --git a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m
index 9f25ebf..ca4185e 100644
--- a/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m	
+++ b/Direct Methods for Solving Linear Systems/Alg4_back_substitution.m	
@@ -1,11 +1,12 @@
 function b = Alg4_back_substitution(U,b)
-% Argorithm 4: Back Substitution (Golub, Loan, Alg. 3.1.2)
-% Returns vetor b with solution to he Ux = b.
+% Argorithm 4: Back Substitution 
+% b = Alg4_back_substitution(U,b) returns vetor b with solution to the
+% Ux = b.
 
 [m, n] = size(U);
 
 if U ~= triu(U)
-    error('Matrix is not upper triangular!')
+    error('Matrix U is not upper triangular!')
 end
 
 if m ~= n
@@ -16,12 +17,14 @@ if length(b) ~= m
     error('Vector b has wrong length!')
 end
 
-% if det(U) < 0.001
-%     error('Matrix is not nonsingular!')
-% end
+if det(U) < eps
+    error('Matrix is not nonsingular!')
+end
+
+
+% The following algorithm is based on the Algrotihm 3.1.2 from [2].
 
 % b(m, :) so that matrices are also accepted
-
 b(m, :) = b(m, :)/U(m, m);
 for i = m-1:-1:1
     b(i, :) = (b(i, :) - U(i, i+1 : m)*b(i+1 : m, :))/U(i, i);
diff --git a/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m b/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m
index f510427..28f598d 100644
--- a/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m	
+++ b/Direct Methods for Solving Linear Systems/Alg5_gauss_jordan_elimination.m	
@@ -1,6 +1,7 @@
 function A = Alg5_gauss_jordan_elimination(A)
 % Algorithm 5: Gauss-Jordan Elimination
-% Argument A is an augmented matrix
+% A = Alg5_gauss_jordan_elimination(A) performs Gauss-Jordan elimination
+% on an augmented matrix A.
 
 [m, n] = size(A);
 
@@ -16,4 +17,5 @@ for k = 1 : m
         end
     end
 end
+
 end
\ No newline at end of file
diff --git a/Direct Methods for Solving Linear Systems/Alg6_RREF.m b/Direct Methods for Solving Linear Systems/Alg6_RREF.m
index 5d055db..c57bb83 100644
--- a/Direct Methods for Solving Linear Systems/Alg6_RREF.m	
+++ b/Direct Methods for Solving Linear Systems/Alg6_RREF.m	
@@ -2,38 +2,38 @@ function A = Alg6_RREF(A)
 % Algorithm 6: Reduced Row Echelon Form (RREF)
 % A = Alg6_RREF(A) returns RREF of matrix A. 
 
-[M, N] = size(A);
+[m, n] = size(A);
 
-n = 0;
+j = 0;
 
-for m = 1 : M
-    n = n + 1;
-    if n > N
+for k = 1 : m
+    j = j + 1;
+    if j > n
         break
     end
     % We want the left-most coefficient to be 1 (pivot)
-    row = A(m, :);
-    if row(m) == 0
-        n = n + 1;
+    row = A(k, :);
+    if row(k) == 0
+        j = j + 1;
     end
-    row = row/row(n);
-    A(m, :) = row;
+    row = row/row(j);
+    A(k, :) = row;
     
-    for i = 1 : M
-        if i ~= m
-            A(i, :) = A(i, :)-(A(i, n))*row;
+    for i = 1 : m
+        if i ~= k
+            A(i, :) = A(i, :)-(A(i, j))*row;
         end
     end
     
-    for i = m + 1 : M
-       A(i:end, m+1:end); % Partial matrix (in which we are looking for non-zero pivots)
-       A(i:end, m+1); % Left-most column
-       if ~any(A(i:end, m+1)) % If the left-most column has only zeros check the next one
-           m = m + 1;
+    for i = k + 1 : m
+       A(i:end, k+1:end); % Partial matrix (in which we are looking for non-zero pivots)
+       A(i:end, k+1); % Left-most column
+       if ~any(A(i:end, k+1)) % If the left-most column has only zeros check the next one
+           k = k + 1;
        end
-       A(i:end, m+1:end);
-       if A(i, m+1) == 0
-           non_zero_row = find(A(i:end,m+1), 1);
+       A(i:end, k+1:end);
+       if A(i, k+1) == 0
+           non_zero_row = find(A(i:end,k+1), 1);
            if isempty(non_zero_row)
                continue
            end
diff --git a/Direct Methods for Solving Linear Systems/README.md b/Direct Methods for Solving Linear Systems/README.md
new file mode 100644
index 0000000..6119207
--- /dev/null
+++ b/Direct Methods for Solving Linear Systems/README.md	
@@ -0,0 +1,8 @@
+# Direct Methods for Solving Linear Systems
+
+## Bibliography
+
+References in the code comments are represented by `[number]`:
+
+  - 1
+  - 2 – Golub, van Loan, Matrix Computations, 3rd edition
\ No newline at end of file