diff --git a/fdm-devito-notebooks/05_nonlin/nonlin_ode.ipynb b/fdm-devito-notebooks/05_nonlin/nonlin_ode.ipynb
index 6a0a6caa..b05bea5c 100644
--- a/fdm-devito-notebooks/05_nonlin/nonlin_ode.ipynb
+++ b/fdm-devito-notebooks/05_nonlin/nonlin_ode.ipynb
@@ -474,11 +474,23 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{dt - \\sqrt{dt^{2} + 4 dt u_{1} - 2 dt + 1} - 1}{2 dt}$"
+      ],
+      "text/plain": [
+       "(dt - sqrt(dt**2 + 4*dt*u_1 - 2*dt + 1) - 1)/(2*dt)"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "import sympy as sym\n",
     "dt, u_1, u = sym.symbols('dt u_1 u')\n",
@@ -488,35 +500,59 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{dt + \\sqrt{dt^{2} + 4 dt u_{1} - 2 dt + 1} - 1}{2 dt}$"
+      ],
+      "text/plain": [
+       "(dt + sqrt(dt**2 + 4*dt*u_1 - 2*dt + 1) - 1)/(2*dt)"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "r2"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-1/dt + 1 - u_1 + dt*(u_1**2 - u_1) + O(dt**2)\n"
+     ]
+    }
+   ],
    "source": [
-    "print r1.series(dt, 0, 2)    # 2 terms in dt, around dt=0"
+    "print(r1.series(dt, 0, 2))    # 2 terms in dt, around dt=0"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "u_1 + dt*(-u_1**2 + u_1) + O(dt**2)\n"
+     ]
+    }
+   ],
    "source": [
-    "print r2.series(dt, 0, 2)"
+    "print(r2.series(dt, 0, 2))"
    ]
   },
   {
@@ -540,24 +576,36 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-oo\n"
+     ]
+    }
+   ],
    "source": [
-    "print r1.limit(dt, 0)"
+    "print(r1.limit(dt, 0))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "u_1\n"
+     ]
+    }
+   ],
    "source": [
-    "print r2.limit(dt, 0)"
+    "print(r2.limit(dt, 0))"
    ]
   },
   {
@@ -624,11 +672,11 @@
    "metadata": {},
    "source": [
     "Since the equation $\\hat F=0$ is only approximate, the solution $u$\n",
-    "does not equal the exact solution $\\uex$ of the exact\n",
-    "equation $F(\\uex)=0$, but we can hope that $u$ is closer to\n",
-    "$\\uex$ than $u^{-}$ is, and hence it makes sense to repeat the\n",
+    "does not equal the exact solution $u_\\text{e}$ of the exact\n",
+    "equation $F(u_\\text{e})=0$, but we can hope that $u$ is closer to\n",
+    "$u_\\text{e}$ than $u^{-}$ is, and hence it makes sense to repeat the\n",
     "procedure, i.e., set $u^{-}=u$ and solve $\\hat F(u)=0$ again.\n",
-    "There is no guarantee that $u$ is closer to $\\uex$ than $u^{-}$,\n",
+    "There is no guarantee that $u$ is closer to $u_\\text{e}$ than $u^{-}$,\n",
     "but this approach has proven to be effective in a wide range of\n",
     "applications.\n",
     "\n",
@@ -784,7 +832,7 @@
     "\n",
     "We shall later refer to the strategy of taking one Picard step, or\n",
     "equivalently, linearizing terms with use of the solution at the\n",
-    "previous time step, as the *Picard1* method. mathcal{I}_t is a widely used\n",
+    "previous time step, as the *Picard1* method. It is a widely used\n",
     "approach in science and technology, but with some limitations if\n",
     "$\\Delta t$ is not sufficiently small (as will be illustrated later).\n",
     "\n",
@@ -950,7 +998,7 @@
     "The geometric mean approximation is often very effective for\n",
     "linearizing quadratic nonlinearities. Both the arithmetic and geometric mean\n",
     "approximations have truncation errors of order $\\Delta t^2$ and are\n",
-    "therefore compatible with the truncation error $\\Oof{\\Delta t^2}$\n",
+    "therefore compatible with the truncation error $\\mathcal{O}(\\Delta t^2)$\n",
     "of the centered difference approximation for $u^\\prime$ in the Crank-Nicolson\n",
     "method.\n",
     "\n",
@@ -1066,7 +1114,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "mathcal{I}_t can be shown that the error in iteration $k+1$ of Newton's method is\n",
+    "It can be shown that the error in iteration $k+1$ of Newton's method is\n",
     "proportional to\n",
     "the square of the error in iteration $k$, a result referred to as\n",
     "*quadratic convergence*. This means that for\n",
@@ -1076,7 +1124,7 @@
     "However, the quadratic convergence appears only if $u^k$ is sufficiently\n",
     "close to the solution. Further away from the solution the method can\n",
     "easily converge very slowly or diverge. The reader is encouraged to do\n",
-    "[nonlin:exer:Newton:problems1](#nonlin:exer:Newton:problems1) to get a better understanding\n",
+    "[Problem 3](nonlin_exer.ipynb#nonlin:exer:Newton:problems1) to get a better understanding\n",
     "for the behavior of the method.\n",
     "\n",
     "Application of Newton's method to the logistic equation discretized\n",
@@ -1229,7 +1277,7 @@
     "<div id=\"nonlin:timediscrete:logistic:impl\"></div>\n",
     "\n",
     "\n",
-    "The program [`logistic.py`](${src_nonlin}/logistic.py) contains\n",
+    "The program [`logistic.py`](https://github.com/devitocodes/devito_book/blob/master/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py) contains\n",
     "implementations of all the methods described above.\n",
     "Below is an extract of the file showing how the Picard and Newton\n",
     "methods are implemented for a Backward Euler discretization of\n",
@@ -1241,9 +1289,7 @@
   {
    "cell_type": "code",
    "execution_count": 7,
-   "metadata": {
-    "collapsed": false
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "def BE_logistic(u0, dt, Nt, choice='Picard',\n",
@@ -1302,9 +1348,7 @@
   {
    "cell_type": "code",
    "execution_count": 8,
-   "metadata": {
-    "collapsed": false
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "def CN_logistic(u0, dt, Nt):\n",
@@ -1336,7 +1380,7 @@
     "Newton's method with one step comes far below the $\\epsilon_r$ tolerance,\n",
     "while the Picard iteration needs on average 7 iterations to bring the\n",
     "residual down to $\\epsilon_r=10^{-1}$, which gives insufficient\n",
-    "accuracy in the solution of the nonlinear equation. mathcal{I}_t is obvious\n",
+    "accuracy in the solution of the nonlinear equation. It is obvious\n",
     "that the Picard1 method gives significant errors in addition to\n",
     "the time discretization unless the time step is as small as in\n",
     "the lower right plot.\n",
@@ -1549,7 +1593,7 @@
     "but with $f(u,t)=\\sin(2(u+1))$, the $f(u^{-},t)u/u^{-}$ trick\n",
     "leads to 7, 9, and 11 iterations during the first three steps, while\n",
     "$f(u^{-},t)$ demands 17, 21, and 20 iterations.\n",
-    "(Experiments can be done with the code [`ODE_Picard_tricks.py`](${src_nonlin}/ODE_Picard_tricks.py).)\n",
+    "(Experiments can be done with the code [`ODE_Picard_tricks.py`](https://github.com/devitocodes/devito_book/blob/master/fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py)\n",
     "\n",
     "\n",
     "\n",
@@ -2327,7 +2371,7 @@
     "can be misleading when the initial start value for the iteration is\n",
     "very close to the solution, since an unnecessary reduction in\n",
     "the error measure is enforced. In such cases the absolute criteria\n",
-    "work better. mathcal{I}_t is common to combine the absolute and relative measures\n",
+    "work better. It is common to combine the absolute and relative measures\n",
     "of the size of the residual,\n",
     "as in"
    ]
@@ -2647,7 +2691,25 @@
    ]
   }
  ],
- "metadata": {},
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.2"
+  }
+ },
  "nbformat": 4,
  "nbformat_minor": 4
 }
diff --git a/fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py b/fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py
index 1ca17e6f..cffbaffd 100644
--- a/fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py
+++ b/fdm-devito-notebooks/05_nonlin/src-nonlin/ODE_Picard_tricks.py
@@ -24,12 +24,12 @@ def f2(u, t):
 tp = linspace(0, 4, 11)
 u1, t = solver1.solve(tp)
 u2, t = solver2.solve(tp)
-print 'Picard it:', solver1.num_iterations_total
-print 'Picard2 it:', solver2.num_iterations_total
+print('Picard it:', solver1.num_iterations_total)
+print('Picard2 it:', solver2.num_iterations_total)
 
 from scitools.std import plot
 plot(t, u1, t, u2, legend=('Picard', 'Picard2'))
-raw_input()
+input()
 
 """
 f(u,t) = -u**3:
diff --git a/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py b/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py
index d5b0f3a6..4ad588b9 100644
--- a/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py
+++ b/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic.py
@@ -1,11 +1,16 @@
 import numpy as np
 
 def FE_logistic(u0, dt, Nt):
-    u = np.zeros(N+1)
-    u[0] = u0
-    for n in range(Nt):
-        u[n+1] = u[n] + dt*(u[n] - u[n]**2)
-    return u
+    t = Dimension('t', spacing=Constant('h_t'))
+    
+    u = TimeFunction(name='u', dimensions=(t,), shape=(Nt+1))
+    u.data[0] = u0
+    
+    pde = u.dtr - u + u**2
+    stencil = Eq(u.forward, solve(pde, u.forward))
+    op = Operator(stencil)
+    op.apply(h_t=dt, t_M=Nt-1)
+    return u.data
 
 def quadratic_roots(a, b, c):
     delta = b**2 - 4*a*c
@@ -77,22 +82,22 @@ def quadratic_root_goes_to_infinity():
         a = dt
         b = 1 - dt
         c = -0.1
-        print dt, quadratic_roots(a, b, c)
+        print(dt, quadratic_roots(a, b, c))
 
 def sympy_analysis():
-    print 'sympy calculations'
+    print('sympy calculations')
     import sympy as sym
     dt, u_1, u = sym.symbols('dt u_1 u')
     r1, r2 = sym.solve(dt*u**2 + (1-dt)*u - u_1, u)
-    print r1
-    print r2
-    print r1.series(dt, 0, 2)
-    print r2.series(dt, 0, 2)
-    print r1.limit(dt, 0)
-    print r2.limit(dt, 0)
+    print(r1)
+    print(r2)
+    print(r1.series(dt, 0, 2))
+    print(r2.series(dt, 0, 2))
+    print(r1.limit(dt, 0))
+    print(r2.limit(dt, 0))
 
 sympy_analysis()
-print '-----------------------------------------------------'
+print('-----------------------------------------------------')
 T = 9
 try:
     dt = float(sys.argv[1])
@@ -113,8 +118,8 @@ def sympy_analysis():
 u_CN = CN_logistic(0.1, dt, N)
 
 from numpy import mean
-print 'Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3)))
-print 'Newton mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE4)))
+print('Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3))))
+print('Newton mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE4))))
 
 t = np.linspace(0, dt*N, N+1)
 plot(t, u_FE, t, u_BE2, t, u_BE3, t, u_BE31, t, u_BE4, t, u_CN,
@@ -125,8 +130,8 @@ def sympy_analysis():
 savefig(filestem + '_u.png')
 savefig(filestem + '_u.pdf')
 figure()
-plot(range(1, len(iter_BE3)+1), iter_BE3, 'r-o',
-     range(1, len(iter_BE4)+1), iter_BE4, 'b-o',
+plot(list(range(1, len(iter_BE3)+1)), iter_BE3, 'r-o',
+     list(range(1, len(iter_BE4)+1)), iter_BE4, 'b-o',
      legend=['Picard', 'Newton'], title='dt=%g, eps=%.0E' % (dt, eps_r),
      axis=[1, N+1, 0, max(iter_BE3 + iter_BE4)+1],
      xlabel='Time level', ylabel='No of iterations')
diff --git a/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic_gen.py b/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic_gen.py
index 3a4864e7..a92dab95 100644
--- a/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic_gen.py
+++ b/fdm-devito-notebooks/05_nonlin/src-nonlin/logistic_gen.py
@@ -77,16 +77,16 @@ def quadratic_root_goes_to_infinity():
         a = dt
         b = 1 - dt
         c = -0.1
-        print dt, quadratic_roots(a, b, c)
+        print(dt, quadratic_roots(a, b, c))
 
-print 'sympy calculations'
+print('sympy calculations')
 import sympy as sym
 dt, u_1, u = sym.symbols('dt u_1 u')
 r1, r2 = sym.solve(dt*u**2 + (1-dt)*u - u_1, u)
-print r1
-print r2
-print r1.series(dt, 0, 2)
-print r2.series(dt, 0, 2)
+print(r1)
+print(r2)
+print(r1.series(dt, 0, 2))
+print(r2.series(dt, 0, 2))
 
 T = 9
 try:
@@ -100,8 +100,8 @@ def quadratic_root_goes_to_infinity():
 N = int(round(T/float(dt)))
 
 u_BE3, iter_BE3 = BE_logistic(0.1, dt, N, 'Picard', eps_r, omega)
-print iter_BE3
-print 'Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3)))
+print(iter_BE3)
+print('Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3))))
 sys.exit(0)
 u_FE = FE_logistic(0.1, dt, N)
 u_BE1, _ = BE_logistic(0.1, dt, N, 'r1')
@@ -112,8 +112,8 @@ def quadratic_root_goes_to_infinity():
 u_CN = CN_logistic(0.1, dt, N)
 
 from numpy import mean
-print 'Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3)))
-print 'Newton mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE4)))
+print('Picard mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE3))))
+print('Newton mean no of iterations (dt=%g):' % dt, int(round(mean(iter_BE4))))
 
 t = np.linspace(0, dt*N, N+1)
 plot(t, u_FE, t, u_BE2, t, u_BE3, t, u_BE31, t, u_BE4, t, u_CN,
@@ -124,8 +124,8 @@ def quadratic_root_goes_to_infinity():
 savefig(filestem + '_u.png')
 savefig(filestem + '_u.pdf')
 figure()
-plot(range(1, len(iter_BE3)+1), iter_BE3, 'r-o',
-     range(1, len(iter_BE4)+1), iter_BE4, 'b-o',
+plot(list(range(1, len(iter_BE3)+1)), iter_BE3, 'r-o',
+     list(range(1, len(iter_BE4)+1)), iter_BE4, 'b-o',
      legend=['Picard', 'Newton'], title='dt=%g, eps=%.0E' % (dt, eps_r),
      axis=[1, N+1, 0, max(iter_BE3 + iter_BE4)+1],
      xlabel='Time level', ylabel='No of iterations')
diff --git a/fdm-devito-notebooks/05_nonlin/src-nonlin/split_diffu_react.py b/fdm-devito-notebooks/05_nonlin/src-nonlin/split_diffu_react.py
index 586051a9..6dda2427 100644
--- a/fdm-devito-notebooks/05_nonlin/src-nonlin/split_diffu_react.py
+++ b/fdm-devito-notebooks/05_nonlin/src-nonlin/split_diffu_react.py
@@ -52,7 +52,7 @@ def diffusion_theta(I, a, f, L, dt, F, t, T, step_no, theta=0.5,
         diagonals=[diagonal, lower, upper],
         offsets=[0, -1, 1], shape=(Nx+1, Nx+1),
         format='csr')
-    #print A.todense()
+    #print(A.todense())
 
     # Allow f to be None or 0
     if f is None or f == 0:
@@ -274,7 +274,7 @@ def action(u, x, t, n):
     for Nx in Nx_values:        
         dx = L/Nx           
         if scheme == 'Strang_splitting_2ndOrder':
-            print 'Strang splitting with 2nd order schemes...'
+            print('Strang splitting with 2nd order schemes...')
             # In this case, E = C*h**r (with r = 2) and since 
             # h = dx = K*dt, the ratio dt/dx must be constant.
             # To fulfill this demand, we must let F change
@@ -283,7 +283,7 @@ def action(u, x, t, n):
             # Initially, we simply choose F = 0.5.
 
             dt = F/a*dx**2
-            #print 'dt/dx:', dt/dx            
+            #print('dt/dx:', dt/dx)
             Nt = int(round(T/float(dt)))
             t = np.linspace(0, Nt*dt, Nt+1)   # global time            
             Strang_splitting_2ndOrder(I=I, a=a, b=b, f=f, L=L, dt=dt, 
@@ -301,39 +301,39 @@ def action(u, x, t, n):
             Nt = int(round(T/float(dt)))
             t = np.linspace(0, Nt*dt, Nt+1)   # global time
             if scheme == 'diffusion':
-                print 'FE on whole eqn...'
+                print('FE on whole eqn...')
                 diffusion_theta(I, a, f, L, dt, F, t, T,
                                 step_no=0, theta=0,
                                 u_L=0, u_R=0, user_action=action)                         
                 h.append(dx**2)
             elif scheme == 'ordinary_splitting':
-                print 'Ordinary splitting...'
+                print('Ordinary splitting...')
                 ordinary_splitting(I=I, a=a, b=b, f=f, L=L, dt=dt, 
                                    dt_Rfactor=1, F=F, t=t, T=T,
                                    user_action=action)        
                 h.append(dx**2)
             elif scheme == 'Strang_splitting_1stOrder':
-                print 'Strang splitting with 1st order schemes...'
+                print('Strang splitting with 1st order schemes...')
                 Strang_splitting_1stOrder(I=I, a=a, b=b, f=f, L=L, dt=dt, 
                                           dt_Rfactor=1, F=F, t=t, T=T,
                                           user_action=action)       
                 h.append(dx**2)
             else:
-                print 'Unknown scheme requested!'
+                print('Unknown scheme requested!')
                 sys.exit(0)
         
-            #print 'dt/dx**2:', dt/dx**2            
+            #print('dt/dx**2:', dt/dx**2)
         
         E.append(error)
         Nx *= 2         # Nx doubled gives dx/2
 
-    print 'E:', E
-    print 'h:', h
+    print('E:', E)
+    print('h:', h)
         
     # Convergence rates 
     r = [np.log(E[i]/E[i-1])/np.log(h[i]/h[i-1])
          for i in range(1,len(Nx_values))]
-    print 'Computed rates:', r
+    print('Computed rates:', r)
 
 
 if __name__ == '__main__':