diff --git a/python/e144.py b/python/e144.py
new file mode 100644
index 0000000..2342a82
--- /dev/null
+++ b/python/e144.py
@@ -0,0 +1,129 @@
+import math
+from dataclasses import dataclass
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+@dataclass
+class P:
+    x: float
+    y: float
+
+
+@dataclass
+class Line:
+    m: float
+    n: float
+
+    def plot(self, x_min: float, x_max: float, color='r'):
+        x = np.linspace(x_min, x_max, 10**5)
+        y = self.m * x + self.n
+        plt.plot(x, y, color, label='') 
+
+
+def line_from_two_points(p1: P, p2: P) -> Line:
+    m = (p2.y - p1.y) / (p2.x - p1.x)
+    n = ((p1.y * p2.x) - (p2.y * p1.x)) / (p2.x - p1.x)
+    return Line(m, n)
+
+
+def line_from_one_point(p: P, m: float) -> Line:
+    n = -m * p.x + p.y
+    return Line(m, n)
+
+
+def reflected_slope(m1, m2):
+    """ Math """
+    alpha1 = math.atan(m1)
+    alpha2 = -math.atan(1/m2)
+    reflected_angle = alpha1 - 2*alpha2
+    m3 = math.tan(reflected_angle)
+    return -m3
+
+
+def get_next_point(start_point: P, line: Line) -> P:
+    # With y = m*x + n into 4x**2 + y**2 = 10 get
+    # quadratic equation x**2 + 2mn * x + (n**2 - 10) = 0.
+
+    # Solve quadratic equation
+    a = 4 + line.m ** 2
+    b = 2 * line.m * line.n
+    c = line.n**2 - 100
+    s = math.sqrt(b**2 - 4*a*c)
+
+    # Pick correct x 
+    EPSILON = 10e-9
+    x1 = (-b + s) / (2*a)
+    x2 = (-b - s) / (2*a)
+    x = x2 if abs(x1 - start_point.x) < EPSILON else x1
+
+    # Calculate corresponding y and return point
+    y = x * line.m + line.n
+    return P(x, y)
+
+
+def calc_slope_ellipse(p: P) -> float:
+    """ Given by Project Euler. """
+    return -4 * p.x / p.y
+
+
+def plot_first_ten():
+    x = np.linspace(-10, 10, 10**5)
+    y_positive = np.sqrt(100 - 4*x**2)
+    y_negative = -np.sqrt(100 - 4*x**2)
+
+    plt.figure(figsize=(6,6)) 
+    plt.plot(x, y_positive, 'b', label='y=sqrt(10-4x^2)') 
+    plt.plot(x, y_negative, 'b', label='y=-sqrt(10-4x^2)') 
+    plt.xlabel('x') 
+    plt.ylabel('y') 
+    plt.title('Plot of the ellipse 4x^2 + y^2 = 10') 
+    plt.axis('equal') 
+
+    current_point = P(0.0, 10.1)
+    next_point = P(1.4, -9.6)
+    line = line_from_two_points(current_point, next_point)
+    next_point_computed = get_next_point(current_point, line)
+    line.plot(0, 1.4)
+    assert(next_point == next_point_computed)
+    current_point = next_point
+
+    for _ in range(10):
+        m1 = line.m
+        m2 = calc_slope_ellipse(current_point)
+        m3 = reflected_slope(m1, m2)
+        line = line_from_one_point(current_point, m3)
+        next_point = get_next_point(current_point, line)
+        x_min = min(current_point.x, next_point.x)
+        x_max = max(current_point.x, next_point.x)
+        line.plot(x_min, x_max, 'g')
+        current_point = next_point
+
+    plt.legend()
+    plt.show()
+
+
+def euler_144():
+    plot_first_ten()
+    current_point = P(0.0, 10.1)
+    next_point = P(1.4, -9.6)
+    line = line_from_two_points(current_point, next_point)
+    current_point = next_point
+    counter = 0
+    while True:
+        counter += 1
+        m1 = line.m
+        m2 = calc_slope_ellipse(current_point)
+        m3 = reflected_slope(m1, m2)
+        line = line_from_one_point(current_point, m3)
+        next_point = get_next_point(current_point, line)
+        current_point = next_point
+        if -0.01 <= current_point.x and current_point.x <= 0.01 and current_point.y > 0:
+            break
+    return counter
+
+
+if __name__ == "__main__":
+    solution = euler_144()
+    print("e144.py: " + str(solution))
+    assert(solution == 354)