diff --git a/CHANGELOG.md b/CHANGELOG.md
index b3a3a23a06..9bcc47938b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -44,6 +44,7 @@
 
 ## Breaking changes
 
+-   Moved `results` to separate repositories ([#761](https://github.com/pybamm-team/PyBaMM/pull/761))
 -   The parameters "Bruggeman coefficient" must now be specified separately as "Bruggeman coefficient (electrolyte)" and "Bruggeman coefficient (electrode)"
 -   The current classes (`GetConstantCurrent`, `GetUserCurrent` and `GetUserData`) have now been removed. Please refer to the [`change-input-current` notebook](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/change-input-current.ipynb) for information on how to specify an input current
 -   Parameter functions must now use pybamm functions instead of numpy functions (e.g. `pybamm.exp` instead of `numpy.exp`), as these are then used to construct the expression tree directly. Generally, pybamm syntax follows numpy syntax; please get in touch if a function you need is missing.
diff --git a/README.md b/README.md
index e302eeb08a..789f4eae29 100644
--- a/README.md
+++ b/README.md
@@ -29,6 +29,8 @@ hosted on [Read The Docs](readthedocs.io). A set of slides giving an overview of
 can be found
 [here](https://github.com/pybamm-team/pybamm_summary_slides/blob/master/pybamm.pdf).
 
+For further examples, see the list of repositories that use PyBaMM [here](https://github.com/pybamm-team/pybamm-example-results)
+
 ## How can I obtain & install PyBaMM?
 
 ### Linux
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000000..c717c5d98a
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,4 @@
+# Examples
+
+A collection of Python scripts and Jupyter notebooks that demonstrate how to use PyBaMM.
+For further examples, see the list of repositories that use PyBaMM [here](https://github.com/pybamm-team/pybamm-example-results)
diff --git a/results/2plus1D/README.md b/results/2plus1D/README.md
deleted file mode 100644
index 1bdf963769..0000000000
--- a/results/2plus1D/README.md
+++ /dev/null
@@ -1,9 +0,0 @@
-# Results: "2plus1D" models
-
-This folder contains the scripts used to generate results for the "2plus1D" papers:
-
-1.  Asymptotic Reduction of a Li-ion Pouch Cell Model: Part 1
-  Scott Marquis, Robert Timms, Valentin Sulzer, Colin Please, S Jon Chapman
-
-2.  Asymptotic Reduction of a Li-ion Pouch Cell Model: Part 2
-  Robert Timms, Scott Marquis, Valentin Sulzer, Colin Please, S Jon Chapman
diff --git a/results/2plus1D/_matplotlibrc b/results/2plus1D/_matplotlibrc
deleted file mode 100644
index 7e8738810d..0000000000
--- a/results/2plus1D/_matplotlibrc
+++ /dev/null
@@ -1,27 +0,0 @@
-# Set a 2/1 aspect ratio with automatic layout.
-figure.figsize: 6.4, 4.
-
-# Use LaTeX fonts
-font.family: serif
-font.serif: ComputerModern
-text.usetex: true
-
-# Backend options (42: use TrueType fonts)
-pdf.fonttype: 42
-ps.fonttype: 42
-
-# Set large font sizes for paper figures
-axes.labelsize: 11
-axes.titlesize: 11
-xtick.labelsize: 11
-ytick.labelsize: 11
-legend.fontsize: 11
-legend.numpoints: 1
-legend.scatterpoints: 1
-
-# Set lines and ticks according to the `seaborn-paper` style sheet
-grid.linewidth: 0.8
-lines.linewidth: 1.4
-patch.linewidth: 0.24
-lines.markersize: 5.6
-lines.markeredgewidth: 0
diff --git a/results/2plus1D/compare_lead_acid_1plus1D.py b/results/2plus1D/compare_lead_acid_1plus1D.py
deleted file mode 100644
index a6b015437e..0000000000
--- a/results/2plus1D/compare_lead_acid_1plus1D.py
+++ /dev/null
@@ -1,73 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level and increase recursion limit
-pybamm.set_logging_level("INFO")
-sys.setrecursionlimit(10000)
-
-# load models
-models = [
-    pybamm.lead_acid.Full(name="1D Full"),
-    pybamm.lead_acid.Composite(name="1D composite"),
-    pybamm.lead_acid.LOQS(name="1D LOQS"),
-    pybamm.lead_acid.Full(
-        {"current collector": "potential pair", "dimensionality": 1}, name="1+1D Full"
-    ),
-    pybamm.lead_acid.Composite(
-        {"current collector": "potential pair", "dimensionality": 1},
-        name="1+1D composite",
-    ),
-    pybamm.lead_acid.LOQS(
-        {"current collector": "potential pair", "dimensionality": 1}, name="1+1D LOQS"
-    ),
-]
-
-# load parameter values and process models
-param = models[0].default_parameter_values
-for model in models:
-    param.process_model(model)
-
-# process geometry and discretise models
-meshes = [None] * len(models)
-for i, model in enumerate(models):
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    var = pybamm.standard_spatial_vars
-    var_pts = {
-        var.x_n: 5,
-        var.x_s: 5,
-        var.x_p: 5,
-        var.r_n: 5,
-        var.r_p: 5,
-        var.y: 5,
-        var.z: 5,
-    }
-    meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve models and process time and voltage for plotting on different meshes
-solutions = [None] * len(models)
-times = [None] * len(models)
-voltages = [None] * len(models)
-t_eval = np.linspace(0, 1, 1000)
-for i, model in enumerate(models):
-    solution = model.default_solver.solve(model, t_eval)
-    solutions[i] = solution
-    times[i] = pybamm.ProcessedVariable(
-        model.variables["Time [h]"], solution.t, solution.y
-    )
-    voltages[i] = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i]
-    )
-
-# plot terminal voltage
-t = np.linspace(0, solution.t[-1], 100)
-for i, model in enumerate(models):
-    plt.plot(times[i](t), voltages[i](t), lw=2, label=model.name)
-plt.xlabel("Time [h]", fontsize=15)
-plt.ylabel("Terminal voltage [V]", fontsize=15)
-plt.legend(fontsize=15)
-plt.show()
diff --git a/results/2plus1D/compare_lead_acid_2plus1D.py b/results/2plus1D/compare_lead_acid_2plus1D.py
deleted file mode 100644
index afd7fdc6de..0000000000
--- a/results/2plus1D/compare_lead_acid_2plus1D.py
+++ /dev/null
@@ -1,75 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level and increase recursion limit
-pybamm.set_logging_level("INFO")
-sys.setrecursionlimit(10000)
-
-# load models
-models = [
-    pybamm.lead_acid.Full(name="1D Full"),
-    pybamm.lead_acid.Composite(name="1D composite"),
-    pybamm.lead_acid.LOQS(name="1D LOQS"),
-    pybamm.lead_acid.Full(
-        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D Full"
-    ),
-    pybamm.lead_acid.Composite(
-        {"current collector": "potential pair", "dimensionality": 2},
-        name="2+1D composite",
-    ),
-    pybamm.lead_acid.LOQS(
-        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D LOQS"
-    ),
-]
-
-# load parameter values and process models
-param = models[0].default_parameter_values
-for model in models:
-    param.process_model(model)
-
-# process geometry and discretise models
-meshes = [None] * len(models)
-for i, model in enumerate(models):
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    var = pybamm.standard_spatial_vars
-    var_pts = {
-        var.x_n: 5,
-        var.x_s: 5,
-        var.x_p: 5,
-        var.r_n: 5,
-        var.r_p: 5,
-        var.y: 5,
-        var.z: 5,
-    }
-    meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve models and process time and voltage for plotting on different meshes
-solutions = [None] * len(models)
-times = [None] * len(models)
-voltages = [None] * len(models)
-t_eval = np.linspace(0, 1, 1000)
-for i, model in enumerate(models):
-    if "2+1D" in model.name:
-        model.use_simplify = False  # simplifying jacobian slow for large systems
-    solution = model.default_solver.solve(model, t_eval)
-    solutions[i] = solution
-    times[i] = pybamm.ProcessedVariable(
-        model.variables["Time [h]"], solution.t, solution.y
-    )
-    voltages[i] = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i]
-    )
-
-# plot terminal voltage
-t = np.linspace(0, solution.t[-1], 100)
-for i, model in enumerate(models):
-    plt.plot(times[i](t), voltages[i](t), lw=2, label=model.name)
-plt.xlabel("Time [h]", fontsize=15)
-plt.ylabel("Terminal voltage [V]", fontsize=15)
-plt.legend(fontsize=15)
-plt.show()
diff --git a/results/2plus1D/compare_lithium_ion_2plus1D.py b/results/2plus1D/compare_lithium_ion_2plus1D.py
deleted file mode 100644
index dfe33bba10..0000000000
--- a/results/2plus1D/compare_lithium_ion_2plus1D.py
+++ /dev/null
@@ -1,96 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level and increase recursion limit
-pybamm.set_logging_level("INFO")
-sys.setrecursionlimit(10000)
-
-# load models
-models = [
-    pybamm.lithium_ion.SPM(name="1D SPM"),
-    pybamm.lithium_ion.SPMe(name="1D SPMe"),
-    pybamm.lithium_ion.DFN(name="1D DFN"),
-    pybamm.lithium_ion.SPM(
-        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPM"
-    ),
-    pybamm.lithium_ion.SPMe(
-        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPMe"
-    ),
-    pybamm.lithium_ion.DFN(
-        {"current collector": "potential pair", "dimensionality": 2}, name="2+1D DFN"
-    ),
-]
-
-# load parameter values
-param = models[0].default_parameter_values
-C_rate = 1
-param.update({"C-rate": C_rate})
-# make current collectors not so conductive, just for illustrative purposes
-param.update(
-    {
-        "Negative current collector conductivity [S.m-1]": 5.96e6,
-        "Positive current collector conductivity [S.m-1]": 3.55e6,
-    }
-)
-
-# process models
-for model in models:
-    param.process_model(model)
-
-# process geometry and discretise models
-meshes = [None] * len(models)
-for i, model in enumerate(models):
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    var = pybamm.standard_spatial_vars
-    var_pts = {
-        var.x_n: 5,
-        var.x_s: 5,
-        var.x_p: 5,
-        var.r_n: 5,
-        var.r_p: 5,
-        var.y: 5,
-        var.z: 5,
-    }
-    meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve models and process time and voltage for plotting on different meshes
-solutions = [None] * len(models)
-times = [None] * len(models)
-voltages = [None] * len(models)
-t_eval = np.linspace(0, 1, 1000)
-for i, model in enumerate(models):
-    if "2+1D" in model.name:
-        model.use_simplify = False  # simplifying jacobian slow for large systems
-    solution = model.default_solver.solve(model, t_eval)
-    solutions[i] = solution
-    times[i] = pybamm.ProcessedVariable(
-        model.variables["Time [h]"], solution.t, solution.y
-    )
-    voltages[i] = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i]
-    )
-
-# plot terminal voltage
-t = np.linspace(0, solution.t[-1], 100)
-for i, model in enumerate(models):
-    plt.plot(times[i](t), voltages[i](t), label=model.name)
-plt.xlabel("Time [h]")
-plt.ylabel("Terminal voltage [V]")
-plt.legend()
-# add C-rate, delta, and alpha to title
-delta = param.evaluate(pybamm.standard_parameters_lithium_ion.delta)
-alpha = param.evaluate(pybamm.standard_parameters_lithium_ion.alpha)
-plt.title(
-    r"C-rate = {:3d}, $\alpha$ = {:.6f} , $\delta$ = {:.6f}".format(
-        C_rate, alpha, delta
-    )
-)
-# save and show
-file_name = "discharge_curve_2plus1D_comparison.eps"
-plt.savefig(file_name, format="eps", dpi=1000)
-plt.show()
diff --git a/results/2plus1D/compare_spmecc.py b/results/2plus1D/compare_spmecc.py
deleted file mode 100644
index 2f4cd2f136..0000000000
--- a/results/2plus1D/compare_spmecc.py
+++ /dev/null
@@ -1,196 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level and increase recursion limit
-pybamm.set_logging_level("INFO")
-sys.setrecursionlimit(10000)
-
-# load current collector and SPMe models
-cc_model = pybamm.current_collector.EffectiveResistance2D()
-spme_av = pybamm.lithium_ion.SPMe(name="Average SPMe")
-spme = pybamm.lithium_ion.SPMe(
-    {"current collector": "potential pair", "dimensionality": 2}, name="2+1D SPMe"
-)
-models = {"Current collector": cc_model, "Average SPMe": spme_av, "2+1D SPMe": spme}
-
-# set parameters based on the spme
-param = spme.default_parameter_values
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {
-    var.x_n: 5,
-    var.x_s: 5,
-    var.x_p: 5,
-    var.r_n: 5,
-    var.r_p: 5,
-    var.y: 5,
-    var.z: 5,
-}
-
-# process model and geometry, and discretise
-meshes = {}
-for name, model in models.items():
-    param.process_model(model)
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    meshes[name] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(meshes[name], model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve models -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solutions = {}
-for name, model in models.items():
-    if name == "Current collector":
-        solutions[name] = model.default_solver.solve(model)
-    else:
-        solutions[name] = model.default_solver.solve(model, t_eval)
-
-# plot terminal voltage
-for name in ["Average SPMe", "2+1D SPMe"]:
-    t, y = solutions[name].t, solutions[name].y
-    model = models[name]
-    time = pybamm.ProcessedVariable(model.variables["Time [h]"], t, y)(t)
-    voltage = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"], t, y, mesh=meshes[name]
-    )(t)
-
-    # add current collector Ohmic losses to average SPMEe to get SPMeCC voltage
-    if model.name == "Average SPMe":
-        current = pybamm.ProcessedVariable(model.variables["Current [A]"], t, y)(t)
-        delta = param.evaluate(pybamm.standard_parameters_lithium_ion.delta)
-        R_cc = param.process_symbol(
-            cc_model.variables["Effective current collector resistance [Ohm]"]
-        ).evaluate(
-            t=solutions["Current collector"].t, y=solutions["Current collector"].y
-        )[
-            0
-        ][
-            0
-        ]
-        cc_ohmic_losses = -delta * current * R_cc
-        voltage = voltage + cc_ohmic_losses
-
-    # plot
-    plt.plot(time, voltage, label=model.name)
-plt.xlabel("Time [h]")
-plt.ylabel("Terminal voltage [V]")
-plt.legend()
-
-
-# plot potentials in current collector
-
-# get processed potentials from SPMeCC
-V_av = pybamm.ProcessedVariable(
-    spme_av.variables["Terminal voltage"],
-    solutions["Average SPMe"].t,
-    solutions["Average SPMe"].y,
-    mesh=meshes["Average SPMe"],
-)
-I_av = pybamm.ProcessedVariable(
-    spme_av.variables["Total current density"],
-    solutions["Average SPMe"].t,
-    solutions["Average SPMe"].y,
-    mesh=meshes["Average SPMe"],
-)
-potentials = cc_model.get_processed_potentials(
-    solutions["Current collector"], meshes["Current collector"], param, V_av, I_av
-)
-phi_s_cn_spmecc = potentials["Negative current collector potential [V]"]
-phi_s_cp_spmecc = potentials["Positive current collector potential [V]"]
-
-# get processed potentials from 2+1D SPMe
-phi_s_cn = pybamm.ProcessedVariable(
-    model.variables["Negative current collector potential [V]"],
-    solutions["2+1D SPMe"].t,
-    solutions["2+1D SPMe"].y,
-    mesh=meshes["2+1D SPMe"],
-)
-phi_s_cp = pybamm.ProcessedVariable(
-    model.variables["Positive current collector potential [V]"],
-    solutions["2+1D SPMe"].t,
-    solutions["2+1D SPMe"].y,
-    mesh=meshes["2+1D SPMe"],
-)
-
-# make plot
-l_y = phi_s_cp.y_sol[-1]
-l_z = phi_s_cp.z_sol[-1]
-y_plot = np.linspace(0, l_y, 21)
-z_plot = np.linspace(0, l_z, 21)
-
-
-def plot(t):
-    plt.subplots(figsize=(15, 8))
-    plt.tight_layout()
-    plt.subplots_adjust(left=-0.1)
-
-    # negative current collector potential
-    plt.subplot(221)
-    phi_s_cn_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=t)),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cn}$")
-    plt.set_cmap("cividis")
-    plt.colorbar(phi_s_cn_plot)
-    plt.subplot(222)
-    phi_s_cn_spmecc_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cn_spmecc(y=y_plot, z=z_plot, t=t)),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cn}$ SPMeCC")
-    plt.set_cmap("cividis")
-    plt.colorbar(phi_s_cn_spmecc_plot)
-
-    # positive current collector potential
-    plt.subplot(223)
-    phi_s_cp_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=t)),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cp}$")
-    plt.set_cmap("viridis")
-    plt.colorbar(phi_s_cp_plot)
-    plt.subplot(224)
-    phi_s_cp_spmecc_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cp_spmecc(y=y_plot, z=z_plot, t=t)),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cp}$ SPMeCC")
-    plt.set_cmap("viridis")
-    plt.colorbar(phi_s_cp_spmecc_plot)
-
-    plt.subplots_adjust(
-        top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5
-    )
-
-
-plot(solutions["2+1D SPMe"].t[-1] / 2)
-plt.show()
diff --git a/results/2plus1D/compare_thermal_lithium_ion_2plus1D.py b/results/2plus1D/compare_thermal_lithium_ion_2plus1D.py
deleted file mode 100644
index c1ac708e98..0000000000
--- a/results/2plus1D/compare_thermal_lithium_ion_2plus1D.py
+++ /dev/null
@@ -1,133 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level and increase recursion limit
-pybamm.set_logging_level("INFO")
-sys.setrecursionlimit(10000)
-
-# load models
-models = [
-    pybamm.lithium_ion.SPM({"thermal": "x-lumped"}, name="1D SPM (lumped)"),
-    pybamm.lithium_ion.SPMe({"thermal": "x-lumped"}, name="1D SPMe (lumped)"),
-    pybamm.lithium_ion.DFN({"thermal": "x-lumped"}, name="1D DFN (lumped)"),
-    pybamm.lithium_ion.SPM({"thermal": "x-full"}, name="1D SPM (full)"),
-    pybamm.lithium_ion.SPMe({"thermal": "x-full"}, name="1D SPMe (full)"),
-    pybamm.lithium_ion.DFN({"thermal": "x-full"}, name="1D DFN (full)"),
-    pybamm.lithium_ion.SPM(
-        {
-            "current collector": "potential pair",
-            "dimensionality": 2,
-            "thermal": "xyz-lumped",
-        },
-        name="2+1D SPM (lumped)",
-    ),
-    pybamm.lithium_ion.SPMe(
-        {
-            "current collector": "potential pair",
-            "dimensionality": 2,
-            "thermal": "xyz-lumped",
-        },
-        name="2+1D SPMe (lumped)",
-    ),
-    pybamm.lithium_ion.DFN(
-        {
-            "current collector": "potential pair",
-            "dimensionality": 2,
-            "thermal": "xyz-lumped",
-        },
-        name="2+1D DFN (lumped)",
-    ),
-    pybamm.lithium_ion.SPM(
-        {
-            "current collector": "potential pair",
-            "dimensionality": 2,
-            "thermal": "x-lumped",
-        },
-        name="2+1D SPM (full)",
-    ),
-    pybamm.lithium_ion.SPMe(
-        {
-            "current collector": "potential pair",
-            "dimensionality": 2,
-            "thermal": "x-lumped",
-        },
-        name="2+1D SPMe (full)",
-    ),
-    pybamm.lithium_ion.DFN(
-        {
-            "current collector": "potential pair",
-            "dimensionality": 2,
-            "thermal": "x-lumped",
-        },
-        name="2+1D DFN (full)",
-    ),
-]
-
-# load parameter values
-param = models[0].default_parameter_values
-
-# process models
-for model in models:
-    param.process_model(model)
-
-# process geometry and discretise models
-meshes = [None] * len(models)
-for i, model in enumerate(models):
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    var = pybamm.standard_spatial_vars
-    var_pts = {
-        var.x_n: 5,
-        var.x_s: 5,
-        var.x_p: 5,
-        var.r_n: 5,
-        var.r_p: 5,
-        var.y: 5,
-        var.z: 5,
-    }
-    meshes[i] = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve models and process time and voltage for plotting on different meshes
-solutions = [None] * len(models)
-times = [None] * len(models)
-voltages = [None] * len(models)
-temperatures = [None] * len(models)
-
-t_eval = np.linspace(0, 1, 1000)
-for i, model in enumerate(models):
-    if "2+1D" in model.name:
-        model.use_simplify = False  # simplifying jacobian slow for large systems
-    solution = model.default_solver.solve(model, t_eval)
-    solutions[i] = solution
-    times[i] = pybamm.ProcessedVariable(
-        model.variables["Time [h]"], solution.t, solution.y
-    )
-    voltages[i] = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=meshes[i]
-    )
-    temperatures[i] = pybamm.ProcessedVariable(
-        model.variables["Volume-averaged cell temperature [K]"],
-        solution.t,
-        solution.y,
-        mesh=meshes[i],
-    )
-
-# plot terminal voltage and temperature
-t = np.linspace(0, solution.t[-1], 100)
-plt.subplot(121)
-for i, model in enumerate(models):
-    plt.plot(times[i](t), voltages[i](t), label=model.name)
-plt.xlabel("Time [h]")
-plt.ylabel("Terminal voltage [V]")
-plt.legend()
-plt.subplot(122)
-for i, model in enumerate(models):
-    plt.plot(times[i](t), temperatures[i](t), label=model.name)
-plt.xlabel("Time [h]")
-plt.ylabel("Temperature [K]")
-plt.tight_layout()
-plt.show()
diff --git a/results/2plus1D/dfn_2plus1D.py b/results/2plus1D/dfn_2plus1D.py
deleted file mode 100644
index 32daf9ed58..0000000000
--- a/results/2plus1D/dfn_2plus1D.py
+++ /dev/null
@@ -1,137 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-# load (2+1D) DFN model
-options = {
-    "current collector": "potential pair",
-    "dimensionality": 2,
-    "thermal": "x-lumped",
-}
-model = pybamm.lithium_ion.DFN(options)
-model.use_simplify = False  # simplifying jacobian slow for large systems
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {
-    var.x_n: 5,
-    var.x_s: 5,
-    var.x_p: 5,
-    var.r_n: 5,
-    var.r_p: 5,
-    var.y: 5,
-    var.z: 5,
-}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = pybamm.IDAKLUSolver().solve(model, t_eval)
-
-# TO DO: 2+1D automated plotting
-phi_s_cn = pybamm.ProcessedVariable(
-    model.variables["Negative current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-phi_s_cp = pybamm.ProcessedVariable(
-    model.variables["Positive current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-T = pybamm.ProcessedVariable(
-    model.variables["X-averaged cell temperature [K]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-l_y = phi_s_cp.y_sol[-1]
-l_z = phi_s_cp.z_sol[-1]
-y_plot = np.linspace(0, l_y, 21)
-z_plot = np.linspace(0, l_z, 21)
-
-
-def plot(t):
-    fig, ax = plt.subplots(figsize=(15, 8))
-    plt.tight_layout()
-    plt.subplots_adjust(left=-0.1)
-
-    # find t index
-    ind = (np.abs(solution.t - t)).argmin()
-
-    # negative current collector potential
-    plt.subplot(131)
-    phi_s_cn_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cn}$ [V]")
-    plt.set_cmap("cividis")
-    plt.colorbar(phi_s_cn_plot)
-
-    # positive current collector potential
-    plt.subplot(132)
-    phi_s_cp_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cp}$ [V]")
-    plt.set_cmap("viridis")
-    plt.colorbar(phi_s_cp_plot)
-
-    # temperature
-    plt.subplot(133)
-    T_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$T$ [K]")
-    plt.set_cmap("inferno")
-    plt.colorbar(T_plot)
-
-    plt.subplots_adjust(
-        top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5
-    )
-    plt.show()
-
-
-plot(solution.t[-1] / 2)
diff --git a/results/2plus1D/set_temperature_spm_1plus1D.py b/results/2plus1D/set_temperature_spm_1plus1D.py
deleted file mode 100644
index c469d76079..0000000000
--- a/results/2plus1D/set_temperature_spm_1plus1D.py
+++ /dev/null
@@ -1,51 +0,0 @@
-#
-# Example of 1+1D SPM where the temperature can be set by the user
-#
-
-import pybamm
-import numpy as np
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-model_options = {
-    "current collector": "potential pair",
-    "dimensionality": 1,
-    "thermal": "x-lumped",
-    "external submodels": ["thermal"],
-}
-model = pybamm.lithium_ion.SPMe(model_options)
-
-var = pybamm.standard_spatial_vars
-z_pts = 20
-var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: z_pts}
-
-sim = pybamm.Simulation(model, var_pts=var_pts, C_rate=2)
-
-# Set the temperature (in dimensionless form)
-# T_av = np.linspace(0, 1, z_pts)[:, np.newaxis]
-
-z = np.linspace(0, 1, z_pts)
-t_eval = np.linspace(0, 0.13, 50)
-# step through the solver, setting the temperature at each timestep
-for i in np.arange(1, len(t_eval) - 1):
-    dt = t_eval[i + 1] - t_eval[i]
-    T_av = (np.sin(2 * np.pi * z) * np.sin(2 * np.pi * 100 * t_eval[i]))[
-        :, np.newaxis
-    ]
-    external_variables = {"X-averaged cell temperature": T_av}
-    sim.step(dt, external_variables=external_variables)
-
-sim.plot(
-    [
-        "Terminal voltage [V]",
-        "X-averaged total heating [W.m-3]",
-        "X-averaged cell temperature [K]",
-        "X-averaged negative particle surface concentration [mol.m-3]",
-        "X-averaged positive particle surface concentration [mol.m-3]",
-        "Negative current collector potential [V]",
-        "Positive current collector potential [V]",
-        "Local voltage [V]",
-    ]
-)
-
diff --git a/results/2plus1D/spm_1plus1D.py b/results/2plus1D/spm_1plus1D.py
deleted file mode 100644
index 993ab68ae9..0000000000
--- a/results/2plus1D/spm_1plus1D.py
+++ /dev/null
@@ -1,63 +0,0 @@
-import pybamm
-import numpy as np
-import sys
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-# load (1+1D) SPMe model
-options = {
-    "current collector": "potential pair",
-    "dimensionality": 1,
-    "thermal": "lumped",
-}
-model = pybamm.lithium_ion.SPM(options)
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-C_rate = 1
-current_1C = 24 * param.process_symbol(pybamm.geometric_parameters.A_cc).evaluate()
-param.update(
-    {
-        "Typical current [A]": C_rate * current_1C,
-        "Initial temperature [K]": 298.15,
-        "Negative current collector conductivity [S.m-1]": 1e7,
-        "Positive current collector conductivity [S.m-1]": 1e7,
-        "Heat transfer coefficient [W.m-2.K-1]": 1,
-    }
-)
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10, var.z: 15}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = model.default_solver.solve(model, t_eval)
-
-# plot
-output_variables = [
-    "X-averaged negative particle surface concentration [mol.m-3]",
-    "X-averaged positive particle surface concentration [mol.m-3]",
-    # "X-averaged cell temperature [K]",
-    "Local potenital difference [V]",
-    "Current collector current density [A.m-2]",
-    "Terminal voltage [V]",
-    "Volume-averaged cell temperature [K]",
-]
-plot = pybamm.QuickPlot(model, mesh, solution, output_variables)
-plot.dynamic_plot()
diff --git a/results/2plus1D/spm_2plus1D.py b/results/2plus1D/spm_2plus1D.py
deleted file mode 100644
index 5f95cd31d0..0000000000
--- a/results/2plus1D/spm_2plus1D.py
+++ /dev/null
@@ -1,137 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level
-pybamm.set_logging_level("DEBUG")
-
-# load (2+1D) SPM model
-options = {
-    "current collector": "potential pair",
-    "dimensionality": 2,
-    "thermal": "x-lumped",
-}
-model = pybamm.lithium_ion.SPM(options)
-model.use_simplify = False  # simplifying jacobian slow for large systems
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {
-    var.x_n: 5,
-    var.x_s: 5,
-    var.x_p: 5,
-    var.r_n: 5,
-    var.r_p: 5,
-    var.y: 5,
-    var.z: 5,
-}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = pybamm.IDAKLUSolver().solve(model, t_eval)
-
-# TO DO: 2+1D automated plotting
-phi_s_cn = pybamm.ProcessedVariable(
-    model.variables["Negative current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-phi_s_cp = pybamm.ProcessedVariable(
-    model.variables["Positive current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-T = pybamm.ProcessedVariable(
-    model.variables["X-averaged cell temperature [K]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-l_y = phi_s_cp.y_sol[-1]
-l_z = phi_s_cp.z_sol[-1]
-y_plot = np.linspace(0, l_y, 21)
-z_plot = np.linspace(0, l_z, 21)
-
-
-def plot(t):
-    fig, ax = plt.subplots(figsize=(15, 8))
-    plt.tight_layout()
-    plt.subplots_adjust(left=-0.1)
-
-    # find t index
-    ind = (np.abs(solution.t - t)).argmin()
-
-    # negative current collector potential
-    plt.subplot(131)
-    phi_s_cn_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cn}$ [V]")
-    plt.set_cmap("cividis")
-    plt.colorbar(phi_s_cn_plot)
-
-    # positive current collector potential
-    plt.subplot(132)
-    phi_s_cp_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cp}$ [V]")
-    plt.set_cmap("viridis")
-    plt.colorbar(phi_s_cp_plot)
-
-    # temperature
-    plt.subplot(133)
-    T_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$T$ [K]")
-    plt.set_cmap("inferno")
-    plt.colorbar(T_plot)
-
-    plt.subplots_adjust(
-        top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5
-    )
-    plt.show()
-
-
-plot(solution.t[-1] / 2)
diff --git a/results/2plus1D/spm_2plus1D_tab_grid.py b/results/2plus1D/spm_2plus1D_tab_grid.py
deleted file mode 100644
index 19835003c5..0000000000
--- a/results/2plus1D/spm_2plus1D_tab_grid.py
+++ /dev/null
@@ -1,138 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-# load (2+1D) SPM model
-options = {
-    "current collector": "potential pair",
-    "dimensionality": 2,
-    "thermal": "x-lumped",
-}
-model = pybamm.lithium_ion.SPM(options)
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {
-    var.x_n: 5,
-    var.x_s: 5,
-    var.x_p: 5,
-    var.r_n: 5,
-    var.r_p: 5,
-    var.y: 5,
-    var.z: 5,
-}
-submesh_types = model.default_submesh_types
-submesh_types["current collector"] = pybamm.ScikitExponential2DSubMesh
-# depnding on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = model.default_solver.solve(model, t_eval)
-
-# TO DO: 2+1D automated plotting
-phi_s_cn = pybamm.ProcessedVariable(
-    model.variables["Negative current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-phi_s_cp = pybamm.ProcessedVariable(
-    model.variables["Positive current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-T = pybamm.ProcessedVariable(
-    model.variables["X-averaged cell temperature [K]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-l_y = phi_s_cp.y_sol[-1]
-l_z = phi_s_cp.z_sol[-1]
-y_plot = np.linspace(0, l_y, 21)
-z_plot = np.linspace(0, l_z, 21)
-
-
-def plot(t):
-    fig, ax = plt.subplots(figsize=(15, 8))
-    plt.tight_layout()
-    plt.subplots_adjust(left=-0.1)
-
-    # find t index
-    ind = (np.abs(solution.t - t)).argmin()
-
-    # negative current collector potential
-    plt.subplot(131)
-    phi_s_cn_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cn}$ [V]")
-    plt.set_cmap("cividis")
-    plt.colorbar(phi_s_cn_plot)
-
-    # positive current collector potential
-    plt.subplot(132)
-    phi_s_cp_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cp}$ [V]")
-    plt.set_cmap("viridis")
-    plt.colorbar(phi_s_cp_plot)
-
-    # temperature
-    plt.subplot(133)
-    T_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$T$ [K]")
-    plt.set_cmap("inferno")
-    plt.colorbar(T_plot)
-
-    plt.subplots_adjust(
-        top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5
-    )
-    plt.show()
-
-
-plot(solution.t[-1] / 2)
diff --git a/results/2plus1D/spme_2plus1D.py b/results/2plus1D/spme_2plus1D.py
deleted file mode 100644
index c7256b993a..0000000000
--- a/results/2plus1D/spme_2plus1D.py
+++ /dev/null
@@ -1,137 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-# load (2+1D) SPMe model
-options = {
-    "current collector": "potential pair",
-    "dimensionality": 2,
-    "thermal": "x-lumped",
-}
-model = pybamm.lithium_ion.SPMe(options)
-model.use_simplify = False  # simplifying jacobian slow for large systems
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {
-    var.x_n: 5,
-    var.x_s: 5,
-    var.x_p: 5,
-    var.r_n: 5,
-    var.r_p: 5,
-    var.y: 5,
-    var.z: 5,
-}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = pybamm.IDAKLUSolver().solve(model, t_eval)
-
-# TO DO: 2+1D automated plotting
-phi_s_cn = pybamm.ProcessedVariable(
-    model.variables["Negative current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-phi_s_cp = pybamm.ProcessedVariable(
-    model.variables["Positive current collector potential [V]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-T = pybamm.ProcessedVariable(
-    model.variables["X-averaged cell temperature [K]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-l_y = phi_s_cp.y_sol[-1]
-l_z = phi_s_cp.z_sol[-1]
-y_plot = np.linspace(0, l_y, 21)
-z_plot = np.linspace(0, l_z, 21)
-
-
-def plot(t):
-    fig, ax = plt.subplots(figsize=(15, 8))
-    plt.tight_layout()
-    plt.subplots_adjust(left=-0.1)
-
-    # find t index
-    ind = (np.abs(solution.t - t)).argmin()
-
-    # negative current collector potential
-    plt.subplot(131)
-    phi_s_cn_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cn(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cn}$ [V]")
-    plt.set_cmap("cividis")
-    plt.colorbar(phi_s_cn_plot)
-
-    # positive current collector potential
-    plt.subplot(132)
-    phi_s_cp_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(phi_s_cp(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$\phi_{s,cp}$ [V]")
-    plt.set_cmap("viridis")
-    plt.colorbar(phi_s_cp_plot)
-
-    # temperature
-    plt.subplot(133)
-    T_plot = plt.pcolormesh(
-        y_plot,
-        z_plot,
-        np.transpose(T(y=y_plot, z=z_plot, t=solution.t[ind])),
-        shading="gouraud",
-    )
-
-    plt.axis([0, l_y, 0, l_z])
-    plt.xlabel(r"$y$")
-    plt.ylabel(r"$z$")
-    plt.title(r"$T$ [K]")
-    plt.set_cmap("inferno")
-    plt.colorbar(T_plot)
-
-    plt.subplots_adjust(
-        top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5
-    )
-    plt.show()
-
-
-plot(solution.t[-1] / 2)
diff --git a/results/2plus1D/spmecc.py b/results/2plus1D/spmecc.py
deleted file mode 100644
index 0939667226..0000000000
--- a/results/2plus1D/spmecc.py
+++ /dev/null
@@ -1,59 +0,0 @@
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-# load current collector and SPMe models
-cell_model = pybamm.lithium_ion.SPMe()
-cc_model = pybamm.current_collector.EffectiveResistance2D()
-models = [cell_model, cc_model]
-
-# set parameters based on the cell model
-param = cell_model.default_parameter_values
-
-# make current collectors not so conductive, just for illustrative purposes
-param.update(
-    {
-        "Negative current collector conductivity [S.m-1]": 5.96e6,
-        "Positive current collector conductivity [S.m-1]": 3.55e6,
-    }
-)
-
-# process model and geometry, and discretise
-for model in models:
-    param.process_model(model)
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-    disc.process_model(model)
-
-
-# solve current collector model
-cc_solution = cc_model.default_solver.solve(cc_model)
-
-# solve SPMe -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = cell_model.default_solver.solve(cell_model, t_eval)
-
-# plot terminal voltage
-t, y = solution.t, solution.y
-time = pybamm.ProcessedVariable(cell_model.variables["Time [h]"], t, y)(t)
-voltage = pybamm.ProcessedVariable(cell_model.variables["Terminal voltage [V]"], t, y)
-current = pybamm.ProcessedVariable(cell_model.variables["Current [A]"], t, y)(t)
-delta = param.evaluate(pybamm.standard_parameters_lithium_ion.delta)
-R_cc = param.process_symbol(
-    cc_model.variables["Effective current collector resistance [Ohm]"]
-).evaluate(t=cc_solution.t, y=cc_solution.y)[0][0]
-cc_ohmic_losses = -delta * current * R_cc
-
-plt.plot(time, voltage(t), label="SPMe")
-plt.plot(time, voltage(t) + cc_ohmic_losses, label="SPMeCC")
-plt.xlabel("Time [h]")
-plt.ylabel("Terminal voltage [V]")
-plt.legend()
-plt.show()
diff --git a/results/2plus1D/user_mesh_spm_1plus1D.py b/results/2plus1D/user_mesh_spm_1plus1D.py
deleted file mode 100644
index 9d479d6a66..0000000000
--- a/results/2plus1D/user_mesh_spm_1plus1D.py
+++ /dev/null
@@ -1,73 +0,0 @@
-import pybamm
-import numpy as np
-import sys
-
-# set logging level
-pybamm.set_logging_level("INFO")
-
-# load (1+1D) SPM model
-options = {
-    "current collector": "potential pair",
-    "dimensionality": 1,
-    "thermal": "lumped",
-}
-model = pybamm.lithium_ion.SPM(options)
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-C_rate = 1
-current_1C = 24 * param.evaluate(pybamm.geometric_parameters.A_cc)
-param.update(
-    {
-        "Typical current [A]": C_rate * current_1C,
-        "Initial temperature [K]": 298.15,
-        "Negative current collector conductivity [S.m-1]": 1e5,
-        "Positive current collector conductivity [S.m-1]": 1e5,
-        "Heat transfer coefficient [W.m-2.K-1]": 1,
-        "Negative tab centre z-coordinate [m]": 0,  # negative tab at bottom
-        "Positive tab centre z-coordinate [m]": 0.137,  # positive tab at top
-    }
-)
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh using user-supplied edges in z
-z_edges = np.array([0, 0.025, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 0.975, 1])
-submesh_types = model.default_submesh_types
-submesh_types["current collector"] = pybamm.MeshGenerator(
-    pybamm.UserSupplied1DSubMesh, submesh_params={"edges": z_edges}
-)
-# Need to make sure var_pts for z is one less than number of edges (variables are
-# evaluated at cell centres)
-npts_z = len(z_edges) - 1
-var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10, var.z: npts_z}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model -- simulate one hour discharge
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval = np.linspace(0, t_end, 120)
-solution = model.default_solver.solve(model, t_eval)
-
-# plot
-output_variables = [
-    "X-averaged negative particle surface concentration [mol.m-3]",
-    "X-averaged positive particle surface concentration [mol.m-3]",
-    # "X-averaged cell temperature [K]",
-    "Local current collector potential difference [V]",
-    "Current collector current density [A.m-2]",
-    "Terminal voltage [V]",
-    "Volume-averaged cell temperature [K]",
-]
-plot = pybamm.QuickPlot(model, mesh, solution, output_variables)
-plot.dynamic_plot()
diff --git a/results/change_settings/change_solver_tolerances.py b/results/change_settings/change_solver_tolerances.py
deleted file mode 100644
index 81261e796d..0000000000
--- a/results/change_settings/change_solver_tolerances.py
+++ /dev/null
@@ -1,52 +0,0 @@
-#
-# Compare solution of li-ion battery models when changing solver tolerances
-#
-import numpy as np
-import pybamm
-
-pybamm.set_logging_level("INFO")
-
-# load model
-model = pybamm.lithium_ion.DFN()
-
-
-# process and discretise
-param = model.default_parameter_values
-param.process_model(model)
-geometry = model.default_geometry
-param.process_geometry(geometry)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# tolerances (rtol, atol)
-tols = [[1e-8, 1e-8], [1e-6, 1e-6], [1e-3, 1e-6], [1e-3, 1e-3]]
-
-# solve model
-solutions = [None] * len(tols)
-voltages = [None] * len(tols)
-voltage_rmse = [None] * len(tols)
-labels = [None] * len(tols)
-t_eval = np.linspace(0, 0.17, 100)
-for i, tol in enumerate(tols):
-    solver = pybamm.ScikitsDaeSolver(rtol=tol[0], atol=tol[1])
-    solutions[i] = solver.solve(model, t_eval)
-    voltages[i] = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"],
-        solutions[i].t,
-        solutions[i].y,
-        mesh=mesh,
-    )(solutions[i].t)
-    voltage_rmse[i] = pybamm.rmse(voltages[0], voltages[i])
-    labels[i] = "rtol = {}, atol = {}".format(tol[0], tol[1])
-
-# print RMSE voltage errors vs tighest tolerance
-for i, tol in enumerate(tols):
-    print(
-        "rtol = {}, atol = {}, solve time = {} s, Voltage RMSE = {}".format(
-            tol[0], tol[1], solutions[i].solve_time, voltage_rmse[i]
-        )
-    )
-# plot
-plot = pybamm.QuickPlot([model] * len(solutions), mesh, solutions, labels=labels)
-plot.dynamic_plot()
diff --git a/results/change_settings/compare_var_pts.py b/results/change_settings/compare_var_pts.py
deleted file mode 100644
index 2153a96f62..0000000000
--- a/results/change_settings/compare_var_pts.py
+++ /dev/null
@@ -1,68 +0,0 @@
-#
-# Compare solution of li-ion battery models when varying the number of grid points
-#
-import numpy as np
-import pybamm
-import matplotlib.pyplot as plt
-
-pybamm.set_logging_level("INFO")
-
-# choose number of points per domain (all domains will have same npts)
-Npts = [30, 20, 10, 5]
-
-# create models
-models = [None] * len(Npts)
-for i, npts in enumerate(Npts):
-    models[i] = pybamm.lithium_ion.DFN(name="Npts = {}".format(npts))
-
-# load parameter values and process models and geometry
-param = models[0].default_parameter_values
-for model in models:
-    param.process_model(model)
-
-# set mesh
-meshes = [None] * len(models)
-
-# create geometry and discretise models
-var = pybamm.standard_spatial_vars
-for i, model in enumerate(models):
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    var_pts = {
-        var.x_n: Npts[i],
-        var.x_s: Npts[i],
-        var.x_p: Npts[i],
-        var.r_n: Npts[i],
-        var.r_p: Npts[i],
-    }
-    meshes[i] = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
-    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
-    disc.process_model(model)
-
-# solve model and plot voltage
-solutions = [None] * len(models)
-voltages = [None] * len(models)
-voltage_rmse = [None] * len(models)
-t_eval = np.linspace(0, 0.17, 100)
-for i, model in enumerate(models):
-    solutions[i] = model.default_solver.solve(model, t_eval)
-    voltages[i] = pybamm.ProcessedVariable(
-        model.variables["Terminal voltage [V]"],
-        solutions[i].t,
-        solutions[i].y,
-        mesh=meshes[i],
-    )(solutions[i].t)
-    voltage_rmse[i] = pybamm.rmse(voltages[0], voltages[i])
-    plt.plot(solutions[i].t, voltages[i], label=model.name)
-
-for i, npts in enumerate(Npts):
-    print(
-        "npts = {}, solve time = {} s, Voltage RMSE = {}".format(
-            npts, solutions[i].solve_time, voltage_rmse[i]
-        )
-    )
-
-plt.xlabel(r"$t$")
-plt.ylabel("Voltage [V]")
-plt.legend()
-plt.show()
diff --git a/results/drive_cycles/US06_simulation.py b/results/drive_cycles/US06_simulation.py
deleted file mode 100644
index 64f47fd5e5..0000000000
--- a/results/drive_cycles/US06_simulation.py
+++ /dev/null
@@ -1,42 +0,0 @@
-#
-# Simulate drive cycle loaded from csv file
-#
-import pybamm
-import numpy as np
-
-# load model
-pybamm.set_logging_level("INFO")
-model = pybamm.lithium_ion.SPMe()
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param["Current function"] = "[current data]US06"
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5}
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# simulate US06 drive cycle
-tau = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_eval = np.linspace(0, 600 / tau, 600)
-
-# need to increase max solver steps if solving DAEs along with an erratic drive cycle
-solver = pybamm.CasadiSolver()
-if isinstance(solver, pybamm.DaeSolver):
-    solver.max_steps = 10000
-
-solution = solver.solve(model, t_eval)
-
-# plot
-plot = pybamm.QuickPlot(model, mesh, solution)
-plot.dynamic_plot()
diff --git a/results/drive_cycles/car_current_simulation.py b/results/drive_cycles/car_current_simulation.py
deleted file mode 100644
index 84d37e4f17..0000000000
--- a/results/drive_cycles/car_current_simulation.py
+++ /dev/null
@@ -1,63 +0,0 @@
-#
-# Simulate user-defined current profile
-#
-import pybamm
-import numpy as np
-
-
-def car_current(t):
-    """
-    Piecewise constant current as a function of time in seconds. This is adapted
-    from the file getCarCurrent.m, which is part of the LIONSIMBA toolbox [1]_.
-
-    References
-    ----------
-    .. [1] M Torchio, L Magni, R Bushan Gopaluni, RD Braatz, and D. Raimondoa.
-           LIONSIMBA: A Matlab framework based on a finite volume model suitable
-           for Li-ion battery design, simulation, and control. Journal of The
-           Electrochemical Society, 163(7):1192-1205, 2016.
-    """
-
-    current = (
-        1 * (t >= 0) * (t <= 50)
-        - 0.5 * (t > 50) * (t <= 60)
-        + 0.5 * (t > 60) * (t <= 210)
-        + 1 * (t > 210) * (t <= 410)
-        + 2 * (t > 410) * (t <= 415)
-        + 1.25 * (t > 415) * (t <= 615)
-        - 0.5 * (t > 615)
-    )
-
-    return current
-
-
-# load model
-pybamm.set_logging_level("INFO")
-model = pybamm.lithium_ion.SPMe()
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param["Current function"] = car_current
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# simulate car current for 30 minutes
-tau = param.process_symbol(
-    pybamm.standard_parameters_lithium_ion.tau_discharge
-).evaluate(0)
-t_eval = np.linspace(0, 1800 / tau, 600)
-solution = model.default_solver.solve(model, t_eval)
-
-# plot
-plot = pybamm.QuickPlot(model, mesh, solution)
-plot.dynamic_plot()
diff --git a/results/drive_cycles/discharge_rest.py b/results/drive_cycles/discharge_rest.py
deleted file mode 100644
index 9089648767..0000000000
--- a/results/drive_cycles/discharge_rest.py
+++ /dev/null
@@ -1,76 +0,0 @@
-#
-# Simulate discharge followed by rest
-#
-import pybamm
-import numpy as np
-import matplotlib.pyplot as plt
-
-# load model
-pybamm.set_logging_level("INFO")
-model = pybamm.lithium_ion.SPMe()
-
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-# solve model during discharge stage (1 hour)
-tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_end = 3600 / tau.evaluate(0)
-t_eval1 = np.linspace(0, t_end, 120)
-solution1 = model.default_solver.solve(model, t_eval1)
-
-# process variables for later plotting
-time1 = pybamm.ProcessedVariable(model.variables["Time [h]"], solution1.t, solution1.y)
-voltage1 = pybamm.ProcessedVariable(
-    model.variables["Terminal voltage [V]"], solution1.t, solution1.y, mesh=mesh
-)
-current1 = pybamm.ProcessedVariable(
-    model.variables["Current [A]"], solution1.t, solution1.y, mesh=mesh
-)
-
-# solve again with zero current, using last step of solution1 as initial conditions
-# update the current to be zero
-param["Current function"] = "[zero]"
-param.update_model(model, disc)
-# Note: need to update model.concatenated_initial_conditions *after* update_model,
-# as update_model updates model.concatenated_initial_conditions, by concatenting
-# the (unmodified) initial conditions for each variable
-model.concatenated_initial_conditions = solution1.y[:, -1][:, np.newaxis]
-
-# simulate 1 hour of rest
-t_start = solution1.t[-1]
-t_end = t_start + 3600 / tau.evaluate(0)
-t_eval2 = np.linspace(t_start, t_end, 120)
-solution2 = model.default_solver.solve(model, t_eval2)
-
-# process variables for later plotting
-time2 = pybamm.ProcessedVariable(model.variables["Time [h]"], solution2.t, solution2.y)
-voltage2 = pybamm.ProcessedVariable(
-    model.variables["Terminal voltage [V]"], solution2.t, solution2.y, mesh=mesh
-)
-current2 = pybamm.ProcessedVariable(
-    model.variables["Current [A]"], solution2.t, solution2.y, mesh=mesh
-)
-
-# plot
-plt.subplot(121)
-plt.plot(time1(t_eval1), voltage1(t_eval1), time2(t_eval2), voltage2(t_eval2))
-plt.xlabel("Time [h]")
-plt.ylabel("Voltage [V]")
-plt.subplot(122)
-z = np.linspace(0, 1, 10)
-plt.plot(time1(t_eval1), current1(t_eval1), time2(t_eval2), current2(t_eval2))
-plt.xlabel("Time [h]")
-plt.ylabel("Current [A]")
-plt.show()
diff --git a/results/drive_cycles/user_sin_current_simulation.py b/results/drive_cycles/user_sin_current_simulation.py
deleted file mode 100644
index c3f56d39b9..0000000000
--- a/results/drive_cycles/user_sin_current_simulation.py
+++ /dev/null
@@ -1,63 +0,0 @@
-#
-# Simulate user-defined current profile which takes parameters
-#
-import pybamm
-import numpy as np
-
-
-# create user-defined function
-def my_fun(t, A, omega):
-    return A * np.sin(2 * np.pi * omega * t)
-
-
-# choose amplitude and frequencies
-A = pybamm.electrical_parameters.I_typ
-frequencies = [0.1, 1]
-
-# load models (need to create new instances of model, not copies)
-pybamm.set_logging_level("INFO")
-models = [None] * len(frequencies)
-for i in range(len(frequencies)):
-    models[i] = pybamm.lithium_ion.SPM()
-
-# load parameter values and process models
-param = models[0].default_parameter_values
-for i, frequency in enumerate(frequencies):
-
-    def current(t):
-        return my_fun(t, A, frequency)
-
-    param.update({"Current function": current})
-    param.process_model(models[i])
-
-# discretise models
-for model in models:
-    # create geometry
-    geometry = model.default_geometry
-    param.process_geometry(geometry)
-    mesh = pybamm.Mesh(
-        geometry, models[-1].default_submesh_types, model.default_var_pts
-    )
-    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-    disc.process_model(model)
-
-# Example: simulate for 30 seconds
-simulation_time = 30  # end time in seconds
-tau = param.process_symbol(
-    pybamm.standard_parameters_lithium_ion.tau_discharge
-).evaluate(0)
-
-# loop over frequencies
-solutions = [None] * len(frequencies)
-labels = [None] * len(frequencies)
-for i, frequency in enumerate(frequencies):
-    # need enough timesteps to resolve output
-    npts = 20 * simulation_time * frequency
-    t_eval = np.linspace(0, simulation_time / tau, npts)
-    solutions[i] = model.default_solver.solve(model, t_eval)
-    labels[i] = "Frequency: {} Hz".format(frequency)
-
-# plot
-output_variables = ["Current [A]", "Terminal voltage [V]"]
-plot = pybamm.QuickPlot(models, mesh, solutions, output_variables, labels)
-plot.dynamic_plot()
diff --git a/tests/unit/test_parameters/test_standard_parameters_lead_acid.py b/tests/unit/test_parameters/test_standard_parameters_lead_acid.py
index 8c7ec90dce..fed248f02f 100644
--- a/tests/unit/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/unit/test_parameters/test_standard_parameters_lead_acid.py
@@ -17,7 +17,7 @@ def test_scipy_constants(self):
     def test_all_defined(self):
         parameters = pybamm.standard_parameters_lead_acid
         parameter_values = pybamm.lead_acid.BaseModel().default_parameter_values
-        output_file = "results/2019_09_sulzer_thesis/parameters.txt"
+        output_file = "lead_acid_parameters.txt"
         pybamm.print_parameters(parameters, parameter_values, output_file)
         # test print_parameters with dict and without C-rate
         del parameter_values["Cell capacity [A.h]"]