diff --git a/covid/python/covid.py b/covid/python/covid.py
new file mode 100644
index 0000000000000000000000000000000000000000..c59dee7c38a7c959009f424c37b96084c184155b
--- /dev/null
+++ b/covid/python/covid.py
@@ -0,0 +1,55 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def s(dt, beta, S0, I0, N):
+    return S0 - dt * (beta * S0 * I0 / N)
+
+def ii(dt, beta, lamb, S0, I0, N):
+    return I0 + dt * (beta * S0 * I0 / N - lamb * I0)
+
+def r(dt, lamb, R0, I0, N):
+    return R0 + dt * lamb * I0
+
+def compute_n(S, R, I):
+    return S + R + I
+
+def timestep(S0, I0, R0, dt, beta, lamb, N):
+    S1 = s(dt, beta, S0, I0, N)
+    I1 = ii(dt, beta, lamb, S0, I0, N)
+    R1 = r(dt, lamb, R0, I0, N)
+    return S1, I1, R1
+
+S0 = 400000
+I0 = 1
+R0 = 0
+max_t = 100
+n_steps = 1000
+dt = max_t / n_steps
+
+
+N = compute_n(S0, I0, R0) 
+
+lamb = 1.0 / 14.0
+beta = 1
+
+s_list = [S0]
+r_list = [R0]
+i_list = [I0]
+t_list = [0]
+for i in range(0, n_steps):
+    S1, I1, R1 = timestep(s_list[i], i_list[i], r_list[i], dt, beta, lamb, N)
+    s_list.append(S1)
+    i_list.append(I1)
+    r_list.append(R1)
+    t_list.append((i+1)*dt)
+
+s = np.array(s_list)
+r = np.array(r_list)
+ii = np.array(i_list)
+t = np.array(t_list)
+
+plt.plot(t, s, 'b')
+plt.plot(t, r, 'r')
+plt.plot(t, ii, 'k')
+plt.show()
+