diff --git a/pySDC/playgrounds/dedalus/README.md b/pySDC/playgrounds/dedalus/README.md
index 2b6ebf4c6..1163e1f5b 100644
--- a/pySDC/playgrounds/dedalus/README.md
+++ b/pySDC/playgrounds/dedalus/README.md
@@ -129,4 +129,6 @@ solver = problem.build_solver(SpectralDeferredCorrectionIMEX)
 timestep = 2e-3  # dummy value for example ...
 while solver.proceed:
     solver.step(timestep)
-```
\ No newline at end of file
+```
+
+A full example script for the 2D Rayleigh-Benard Convection problem can be found [here](./rayleighBenardSDC.py).
\ No newline at end of file
diff --git a/pySDC/playgrounds/dedalus/rayleighBenardSDC.py b/pySDC/playgrounds/dedalus/rayleighBenardSDC.py
new file mode 100644
index 000000000..ca09c5e81
--- /dev/null
+++ b/pySDC/playgrounds/dedalus/rayleighBenardSDC.py
@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Dedalus script simulating 2D horizontally-periodic Rayleigh-Benard convection
+using Spectral Deferred Corrections.
+This script demonstrates solving a 2D Cartesian initial value problem. It can
+be ran serially or in parallel, and uses the built-in analysis framework to save
+data snapshots to HDF5 files. The `plot_snapshots.py` script can be used to
+produce plots from the saved data. It should take about 5 cpu-minutes to run.
+
+The problem is non-dimensionalized using the box height and freefall time, so
+the resulting thermal diffusivity and viscosity are related to the Prandtl
+and Rayleigh numbers as:
+
+    kappa = (Rayleigh * Prandtl)**(-1/2)
+    nu = (Rayleigh / Prandtl)**(-1/2)
+
+For incompressible hydro with two boundaries, we need two tau terms for each the
+velocity and buoyancy. Here we choose to use a first-order formulation, putting
+one tau term each on auxiliary first-order gradient variables and the others in
+the PDE, and lifting them all to the first derivative basis. This formulation puts
+a tau term in the divergence constraint, as required for this geometry.
+
+To run using e.g. 4 processes:
+    $ mpiexec -n 4 python rayleighBenardSDC.py
+
+Then use https://github.com/DedalusProject/dedalus/blob/master/examples/ivp_2d_rayleigh_benard/plot_snapshots.py
+for plotting.
+"""
+import numpy as np
+import dedalus.public as d3
+from pySDC.playgrounds.dedalus.sdc import SpectralDeferredCorrectionIMEX
+import logging
+logger = logging.getLogger(__name__)
+
+useSDC = True
+SpectralDeferredCorrectionIMEX.setParameters(
+    nSweeps=4,
+    nNodes=4,
+    implSweep="MIN-SR-FLEX",
+    explSweep="PIC")
+
+# Parameters
+Lx, Lz = 4, 1
+Nx, Nz = 512, 128
+Rayleigh = 2e6
+Prandtl = 1
+dealias = 3/2
+stop_sim_time = 50
+timestepper = SpectralDeferredCorrectionIMEX if useSDC else d3.RK222
+timestep = 1e-2/4
+dtype = np.float64
+
+# Bases
+coords = d3.CartesianCoordinates('x', 'z')
+dist = d3.Distributor(coords, dtype=dtype)
+xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias)
+zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz), dealias=dealias)
+
+# Fields
+p = dist.Field(name='p', bases=(xbasis,zbasis))
+b = dist.Field(name='b', bases=(xbasis,zbasis))
+u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))
+tau_p = dist.Field(name='tau_p')
+tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
+tau_b2 = dist.Field(name='tau_b2', bases=xbasis)
+tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
+tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)
+
+# Substitutions
+kappa = (Rayleigh * Prandtl)**(-1/2)
+nu = (Rayleigh / Prandtl)**(-1/2)
+x, z = dist.local_grids(xbasis, zbasis)
+ex, ez = coords.unit_vector_fields(dist)
+lift_basis = zbasis.derivative_basis(1)
+lift = lambda A: d3.Lift(A, lift_basis, -1)
+grad_u = d3.grad(u) + ez*lift(tau_u1) # First-order reduction
+grad_b = d3.grad(b) + ez*lift(tau_b1) # First-order reduction
+
+# Problem
+# First-order form: "div(f)" becomes "trace(grad_f)"
+# First-order form: "lap(f)" becomes "div(grad_f)"
+problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
+problem.add_equation("trace(grad_u) + tau_p = 0")
+problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)")
+problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - u@grad(u)")
+problem.add_equation("b(z=0) = Lz")
+problem.add_equation("u(z=0) = 0")
+problem.add_equation("b(z=Lz) = 0")
+problem.add_equation("u(z=Lz) = 0")
+problem.add_equation("integ(p) = 0") # Pressure gauge
+
+# Solver
+solver = problem.build_solver(timestepper)
+solver.stop_sim_time = stop_sim_time
+
+# Initial conditions
+b.fill_random('g', seed=42, distribution='normal', scale=1e-3) # Random noise
+b['g'] *= z * (Lz - z) # Damp noise at walls
+b['g'] += Lz - z # Add linear background
+
+# Analysis
+snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50)
+snapshots.add_task(b, name='buoyancy')
+snapshots.add_task(-d3.div(d3.skew(u)), name='vorticity')
+
+# Main loop
+try:
+    logger.info('Starting main loop')
+    while solver.proceed:
+        solver.step(timestep)
+        if (solver.iteration-1) % 10 == 0:
+            logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep))
+except:
+    logger.error('Exception raised, triggering end of main loop.')
+    raise
+finally:
+    solver.log_stats()