diff --git a/SteinsGate2D.py b/SteinsGate2D.py
new file mode 100644
index 0000000..d034856
--- /dev/null
+++ b/SteinsGate2D.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+SteinsGate2D.py
+Created on Sat March 6, 2021
+
+@author: David Nolte
+Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)
+
+2D simulation of Stein's Gate Divergence Meter
+"""
+import numpy as np
+from scipy import integrate
+from matplotlib import pyplot as plt
+
+plt.close('all')
+
+def solve_flow(param,lim = [-6,6,-6,6],max_time=20.0):
+
+
+    def flow_deriv(x_y, t0, alpha, beta, gamma):
+        #"""Compute the time-derivative ."""
+        x, y = x_y
+        
+        w = 1
+        R2 = x**2 + y**2
+        R = np.sqrt(R2)
+        arg = (R-2)/0.1
+        env1 = 1/(1+np.exp(arg))
+        env2 = 1 - env1
+        
+        f = env2*(x*(1/(R-1.99)**2 + 1e-2) - x) + env1*(w*y + w*x*(1 - R))
+        g = env2*(y*(1/(R-1.99)**2 + 1e-2) + y) + env1*(-w*x + w*y*(1 - R))
+        
+        return [f,g]
+    model_title = 'Steins Gate'
+
+    plt.figure()
+    xmin = lim[0]
+    xmax = lim[1]
+    ymin = lim[2]
+    ymax = lim[3]
+    plt.axis([xmin, xmax, ymin, ymax])
+
+    N = 24*4 + 47
+    x0 = np.zeros(shape=(N,2))
+    ind = -1
+    for i in range(0,24):
+        ind = ind + 1
+        x0[ind,0] = xmin + (xmax-xmin)*i/23
+        x0[ind,1] = ymin
+        ind = ind + 1
+        x0[ind,0] = xmin + (xmax-xmin)*i/23
+        x0[ind,1] = ymax
+        ind = ind + 1
+        x0[ind,0] = xmin
+        x0[ind,1] = ymin + (ymax-ymin)*i/23
+        ind = ind + 1
+        x0[ind,0] = xmax
+        x0[ind,1] = ymin + (ymax-ymin)*i/23
+    ind = ind + 1
+    x0[ind,0] = 0.05
+    x0[ind,1] = 0.05
+    
+    for thetloop in range(0,10):
+        ind = ind + 1
+        theta = 2*np.pi*(thetloop)/10
+        ys = 0.125*np.sin(theta)
+        xs = 0.125*np.cos(theta)
+        x0[ind,0] = xs
+        x0[ind,1] = ys
+
+    for thetloop in range(0,10):
+        ind = ind + 1
+        theta = 2*np.pi*(thetloop)/10
+        ys = 1.7*np.sin(theta)
+        xs = 1.7*np.cos(theta)
+        x0[ind,0] = xs
+        x0[ind,1] = ys
+
+    for thetloop in range(0,20):
+        ind = ind + 1
+        theta = 2*np.pi*(thetloop)/20
+        ys = 2*np.sin(theta)
+        xs = 2*np.cos(theta)
+        x0[ind,0] = xs
+        x0[ind,1] = ys
+        
+    ind = ind + 1
+    x0[ind,0] = -3
+    x0[ind,1] = 0.05
+    ind = ind + 1
+    x0[ind,0] = -3
+    x0[ind,1] = -0.05
+    ind = ind + 1
+    x0[ind,0] = 3
+    x0[ind,1] = 0.05
+    ind = ind + 1
+    x0[ind,0] = 3
+    x0[ind,1] = -0.05
+    ind = ind + 1
+    x0[ind,0] = -6
+    x0[ind,1] = 0.00
+    ind = ind + 1
+    x0[ind,0] = 6
+    x0[ind,1] = 0.00
+           
+    colors = plt.cm.prism(np.linspace(0, 1, N))
+                        
+    # Solve for the trajectories
+    t = np.linspace(0, max_time, int(250*max_time))
+    x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t, param)
+                      for x0i in x0])
+
+    for i in range(N):
+        x, y = x_t[i,:,:].T
+        lines = plt.plot(x, y, '-', c=colors[i])
+        plt.setp(lines, linewidth=1)
+
+        
+    plt.show()
+    plt.title(model_title)
+    
+    
+    return t, x_t
+
+param = (0.02,0.5,0.2)        # Steins Gate
+lim = (-6,6,-6,6)
+
+t, x_t = solve_flow(param,lim)
+
+plt.savefig('Steins Gate')