diff --git a/main.py b/main.py
index 220e493..2f70c0b 100644
--- a/main.py
+++ b/main.py
@@ -35,15 +35,51 @@ def interior_point(
         return Result(State.UNAPPLICABLE)
     if (not maximizing):
         C = -C
-    x = x_0
-    solved = False
-    while not solved:
-        # TODO Algorithm steps (refer to numpy.linalg for matrix stuff)
-        pass
+        
+    m = len(A)
+    n = len(A[0])
 
-    # return value (include check for minimization)
-    return Result(State.SOLVED, ...)
+    x = np.ones(n)
+    s = np.ones(m)
 
+    iteration = 0
+
+    while(True):
+
+        for i in range(m):
+            slack = b[i]
+            for j in range(n):
+                slack -= A[i][j] * x[j]
+            s[i] = slack
+        
+        for i in range(min(m, n)):
+            x[i] = s[i]
+
+        D = np.diag(s)
+
+        x_star = np.dot(np.linalg.inv(D), x)
+        A_star = np.dot(A, D)
+        C_star = np.dot(D, C)
+
+        I = np.eye(n)
+        A_star_transpose = np.transpose(A_star)
+        P = I - np.dot(A_star_transpose, np.linalg.inv(np.dot(A_star, A_star_transpose)))
+        P = np.dot(P, A_star)
+
+        C_p = np.dot(P, C_star)
+        Mu = np.max(np.absolute(C_p))
+
+        if Mu < eps:
+            result = np.dot(C, x)
+            return Result(State.SOLVED, objective_function_value=result, solution=x)
+
+        iteration += 1
+
+        if iteration >= 1000:
+            return Result(State.UNSOLVED)
+
+        x_star += (alpha / Mu) * C_p
+        x = np.dot(D, x_star)
 
 # TODO 5 tests (from assignment 1) and comparison with simplex and alpha = 0.9
 def TEST_CASE_GENERAL():