diff --git a/Kuramoto.py b/Kuramoto.py
new file mode 100644
index 0000000..e0f931d
--- /dev/null
+++ b/Kuramoto.py
@@ -0,0 +1,154 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Sat May 11 08:56:41 2019
+
+@author: nolte
+
+D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019)
+"""
+
+# https://www.python-course.eu/networkx.php
+# https://networkx.github.io/documentation/stable/tutorial.html
+# https://networkx.github.io/documentation/stable/reference/functions.html
+
+import numpy as np
+from scipy import integrate
+from matplotlib import pyplot as plt
+import networkx as nx
+from UserFunction import linfit
+import time
+
+tstart = time.time()
+
+plt.close('all')
+
+Nfac = 25   # 25
+N = 50      # 50
+width = 0.2
+
+# function: omegout, yout = coupleN(G)
+def coupleN(G):
+
+    # function: yd = flow_deriv(x_y)
+    def flow_deriv(y,t0):
+                
+        yp = np.zeros(shape=(N,))
+        for omloop  in range(N):
+            temp = omega[omloop]
+            linksz = G.node[omloop]['numlink']
+            for cloop in range(linksz):
+                cindex = G.node[omloop]['link'][cloop]
+                g = G.node[omloop]['coupling'][cloop]
+
+                temp = temp + g*np.sin(y[cindex]-y[omloop])
+            
+            yp[omloop] = temp
+        
+        yd = np.zeros(shape=(N,))
+        for omloop in range(N):
+            yd[omloop] = yp[omloop]
+        
+        return yd
+    # end of function flow_deriv(x_y)
+
+    mnomega = 1.0
+    
+    for nodeloop in range(N):
+        omega[nodeloop] = G.node[nodeloop]['element']
+    
+    x_y_z = omega    
+    
+    # Settle-down Solve for the trajectories
+    tsettle = 100
+    t = np.linspace(0, tsettle, tsettle)
+    x_t = integrate.odeint(flow_deriv, x_y_z, t)
+    x0 = x_t[tsettle-1,0:N]
+    
+    t = np.linspace(1,1000,1000)
+    y = integrate.odeint(flow_deriv, x0, t)
+    siztmp = np.shape(y)
+    sy = siztmp[0]
+        
+    # Fit the frequency
+    m = np.zeros(shape = (N,))
+    w = np.zeros(shape = (N,))
+    mtmp = np.zeros(shape=(4,))
+    btmp = np.zeros(shape=(4,))
+    for omloop in range(N):
+        
+        if np.remainder(sy,4) == 0:
+            mtmp[0],btmp[0] = linfit(t[0:sy//2],y[0:sy//2,omloop]);
+            mtmp[1],btmp[1] = linfit(t[sy//2+1:sy],y[sy//2+1:sy,omloop]);
+            mtmp[2],btmp[2] = linfit(t[sy//4+1:3*sy//4],y[sy//4+1:3*sy//4,omloop]);
+            mtmp[3],btmp[3] = linfit(t,y[:,omloop]);
+        else:
+            sytmp = 4*np.floor(sy/4);
+            mtmp[0],btmp[0] = linfit(t[0:sytmp//2],y[0:sytmp//2,omloop]);
+            mtmp[1],btmp[1] = linfit(t[sytmp//2+1:sytmp],y[sytmp//2+1:sytmp,omloop]);
+            mtmp[2],btmp[2] = linfit(t[sytmp//4+1:3*sytmp/4],y[sytmp//4+1:3*sytmp//4,omloop]);
+            mtmp[3],btmp[3] = linfit(t[0:sytmp],y[0:sytmp,omloop]);
+
+        
+        #m[omloop] = np.median(mtmp)
+        m[omloop] = np.mean(mtmp)
+        
+        w[omloop] = mnomega + m[omloop]
+     
+    omegout = m
+    yout = y
+    
+    return omegout, yout
+    # end of function: omegout, yout = coupleN(G)
+
+
+
+Nlink = N*(N-1)//2      
+omega = np.zeros(shape=(N,))
+omegatemp = width*(np.random.rand(N)-1)
+meanomega = np.mean(omegatemp)
+omega = omegatemp - meanomega
+sto = np.std(omega)
+
+#nodecouple = nx.complete_graph(N)
+nodecouple = nx.erdos_renyi_graph(N,0.1)
+
+lnk = np.zeros(shape = (N,), dtype=int)
+for loop in range(N):
+    nodecouple.node[loop]['element'] = omega[loop]
+    nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop))
+    nodecouple.node[loop]['numlink'] = np.size(list(nx.neighbors(nodecouple,loop)))
+    lnk[loop] = np.size(list(nx.neighbors(nodecouple,loop)))
+
+avgdegree = np.mean(lnk)
+mnomega = 1
+
+facval = np.zeros(shape = (Nfac,))
+yy = np.zeros(shape=(Nfac,N))
+xx = np.zeros(shape=(Nfac,))
+for facloop in range(Nfac):
+    print(facloop)
+    facoef = 0.2
+
+    fac = facoef*(16*facloop/(Nfac))*(1/(N-1))*sto/mnomega
+    for nodeloop in range(N):
+        nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[nodeloop],))
+        for linkloop in range (lnk[nodeloop]):
+            nodecouple.node[nodeloop]['coupling'][linkloop] = fac
+
+    facval[facloop] = fac*avgdegree
+    
+    omegout, yout = coupleN(nodecouple)                           # Here is the subfunction call for the flow
+
+    for omloop in range(N):
+        yy[facloop,omloop] = omegout[omloop]
+
+    xx[facloop] = facval[facloop]
+
+plt.figure(1)
+lines = plt.plot(xx,yy)
+plt.setp(lines, linewidth=0.5)
+plt.show()
+
+elapsed_time = time.time() - tstart
+print('elapsed time = ',format(elapsed_time,'.2f'),'secs')