diff --git a/scripts/adv_disp.py b/scripts/adv_disp.py
index 1ec936a..86ebe84 100644
--- a/scripts/adv_disp.py
+++ b/scripts/adv_disp.py
@@ -1,11 +1,11 @@
 import numpy as np
-from scipy.sparse import linalg
 from pylab import rcParams
 import matplotlib.pyplot as plt
 from matplotlib.patches import Polygon
 from subprocess import call
 import sympy
 import warnings
+from scipy.sparse import linalg
 
 def findomega(Z):
   omega = sympy.Symbol('omega')
diff --git a/scripts/plot_dispersion.py b/scripts/plot_dispersion.py
index 245d7d3..98e5ad1 100644
--- a/scripts/plot_dispersion.py
+++ b/scripts/plot_dispersion.py
@@ -65,8 +65,7 @@ def normalise(R, T, target):
       u0      = solution_linear(u0_val, np.array([[symb]],dtype='complex'))
       ucoarse = solution_linear(u0_val, np.array([[symb_coarse]],dtype='complex'))
       para = parareal(0.0, Tend, nslices, intexact, impeuler, nfine, ncoarse, 0.0, niter_v[0], u0)
-      
-      
+            
       # get update matrix for imp Euler over one slice
       stab_fine   = para.timemesh.slices[0].get_fine_update_matrix(u0)
       stab_fine   = stab_fine
diff --git a/src/trapezoidal.py b/src/trapezoidal.py
new file mode 100644
index 0000000..3c4b363
--- /dev/null
+++ b/src/trapezoidal.py
@@ -0,0 +1,41 @@
+from integrator import integrator
+from solution import solution
+import copy
+
+# Used to compute update matrices
+from solution_linear import solution_linear
+from scipy import sparse
+from scipy import linalg
+import numpy as np
+
+class trapezoidal(integrator):
+
+  def __init__(self, tstart, tend, nsteps):
+    super(trapezoidal, self).__init__(tstart, tend, nsteps)
+    self.order = 2
+
+  def run(self, u0):
+    assert isinstance(u0, solution), "Initial value u0 must be an object of type solution"
+    for i in range(0,self.nsteps):
+      utemp = copy.deepcopy(u0)
+      utemp.f()
+      u0.applyM()
+      u0.axpy(0.5*self.dt, utemp)  
+      u0.solve(0.5*self.dt)
+
+  #
+  # For linear problems My' = A', a backward Euler update corresponds
+  # to 
+  # u_n+1 = (M - dt*A)^(-1)*M*u_n
+  # so that a full update is [(M-dt*A)^(-1)*M]^nsteps
+  def get_update_matrix(self, sol):
+    assert isinstance(sol, solution_linear), "Update function can only be computed for solutions of type solution_linear"
+    # Fetch matrices from solution and make sure they are sparse
+    M = sparse.csc_matrix(sol.M)
+    A = sparse.csc_matrix(sol.A)
+    Rmat = sparse.linalg.inv(M - 0.5*self.dt*A)
+    # this is call is necessary because if Rmat has only 1 entry, it gets converted to a dense array here
+    Rmat = sparse.csc_matrix(Rmat)
+    Rmat = Rmat.dot(M+0.5*self.dt*A)
+    Rmat = Rmat**self.nsteps
+    return Rmat
diff --git a/tests/test_trapezoidal.py b/tests/test_trapezoidal.py
new file mode 100644
index 0000000..9d798a5
--- /dev/null
+++ b/tests/test_trapezoidal.py
@@ -0,0 +1,60 @@
+import sys
+sys.path.append('../src')
+
+import numpy as np
+from scipy import sparse
+from scipy import linalg
+from integrator import integrator
+from trapezoidal import trapezoidal
+from solution_linear import solution_linear
+from nose.tools import *
+import unittest
+
+class TestTrapezoidal(unittest.TestCase):
+  
+  def setUp(self):
+    self.ndof = np.random.randint(255)
+    self.A    = sparse.spdiags([ np.ones(self.ndof), -2.0*np.ones(self.ndof), np.ones(self.ndof)], [-1,0,1], self.ndof, self.ndof, format="csc")
+    self.M    = sparse.spdiags([ np.random.rand(self.ndof) ], [0], self.ndof, self.ndof, format="csc")
+    self.sol  = solution_linear(np.ones(self.ndof), self.A, self.M)
+
+  # Can instantiate object
+  def test_caninstantiate(self):
+    tp = trapezoidal(0.0, 1.0, 10)
+
+  # Throws if tend < tstart
+  def test_tendbeforetstartthrows(self):
+    with self.assertRaises(AssertionError):
+      tp = trapezoidal(1.0, 0.0, 10)
+
+  # See if run function can be called
+  def test_cancallrun(self):
+    tp = trapezoidal(0.0, 1.0, 10)
+    tp.run(self.sol)
+
+  # See if run does the same as the update matrix for a scalar problem
+  def test_callcorrectscalar(self):
+    eig = -1.0
+    u0 = solution_linear(np.array([1.0]), np.array([[eig]]))
+    nsteps = 50
+    tp = trapezoidal(0.0, 1.0, nsteps)
+    tp.run(u0)
+    assert abs(u0.y - np.exp(-1.0))<5e-3, ("Very wrong solution. Error: %5.2e" % abs(u0.y - np.exp(-1.0)))
+    Rmat = tp.get_update_matrix(u0)
+    Rmat_ie = Rmat[0,0]
+    Rmat_ex = ((1.0 + 0.5*tp.dt*eig)/(1.0 - 0.5*tp.dt*eig))**nsteps
+    assert abs(Rmat_ie - Rmat_ex)<1e-14, ("Update function generated by trapezoidal rule for scalar case does not match exact value. Error: %5.3e" % abs(Rmat_ie - Rmat_ex))
+
+  # See if run does the same as the update matrix
+  def test_callcorrect(self):
+    u0 = solution_linear(np.ones(self.ndof), self.A, self.M)
+    nsteps = 13
+    tp = trapezoidal(0.0, 1.0, nsteps)
+    tp.run(u0)
+    # Compute output through update matrix and compare
+    Rmat = tp.get_update_matrix(u0)
+    yend = Rmat.dot(np.ones((self.ndof,1)))
+    sol_end = solution_linear( yend, self.A, self.M )
+    sol_end.axpy(-1.0, u0)
+    assert sol_end.norm()<2e-12, ("Output from trapezoidal rule integrator differs from result computed with power of update matrix -- norm of difference: %5.3e" % sol_end.norm())
+