diff --git a/Direct Methods for Solving Linear Systems/back_substitution.m b/Direct Methods for Solving Linear Systems/back_substitution.m
new file mode 100644
index 0000000..2d1a1f0
--- /dev/null
+++ b/Direct Methods for Solving Linear Systems/back_substitution.m	
@@ -0,0 +1,21 @@
+function back_substitution(U,b)
+
+[n, m] = size(U);
+if n ~= m
+    error('Matrix is not squared!')
+end
+
+if length(b) ~= n
+    error('Vector b has wrong length!')
+end
+
+if det(A) == 0
+    error('Matrix is not nonsingular!') 
+end
+
+b(n) = b(n)/U(n, n)
+for i = n-1:-1:n
+    
+end
+
+end
\ No newline at end of file
diff --git a/Direct Methods for Solving Linear Systems/main.m b/Direct Methods for Solving Linear Systems/main.m
new file mode 100644
index 0000000..37da31a
--- /dev/null
+++ b/Direct Methods for Solving Linear Systems/main.m	
@@ -0,0 +1,8 @@
+clear all;
+
+B = [2, -1, 0, 0; -1, 2, -1, 0;0, -1, 2, -1; 0, 0, -1, 2]
+
+b = [0;0;0;5]
+
+[U, Mk] = outer_product_gaussian_elimination(B)
+back_substitution(U, b)
\ No newline at end of file
diff --git a/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m b/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m
index 012c1de..740deec 100644
--- a/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m	
+++ b/Direct Methods for Solving Linear Systems/outer_product_gaussian_elimination.m	
@@ -1,5 +1,5 @@
 % Outer Product Gaussian Elimination (Alg. 3.2.1)
-function outer_product_gaussian_elimination(A)
+function U, Mk = outer_product_gaussian_elimination(A)
 
 [n, m] = size(A);
 if n ~= m
@@ -13,8 +13,11 @@ end
  A
  
  for k = 1 : n-1
-     k
-     rows = k + 1 : n
-     A(rows, k) = A(rows, k)/A(k, k)
-     A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows)
- end
\ No newline at end of file
+     rows = k + 1 : n;
+     A(rows, k) = A(rows, k)/A(k, k);
+     A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
+     A
+ end
+ 
+U = triu(A)
+Mk = diag(A,-1) % Gauss vector