diff --git a/.gitignore b/.gitignore
index 4c691e285..e567fff73 100644
--- a/.gitignore
+++ b/.gitignore
@@ -45,6 +45,7 @@ dist/
 # coverage
 .coverage
 coverage.xml
+htmlcov/
 
 # black setup file seems to make Travis CI fail
 pyproject.toml
@@ -61,6 +62,3 @@ pyvenv.cfg
 
 # sundials
 sundials
-
-# code coverage
-htmlcov/
diff --git a/examples/scripts/compare_comsol/compare_comsol_DFN.py b/examples/scripts/compare_comsol/compare_comsol_DFN.py
index ba34e461e..0a6b982c9 100644
--- a/examples/scripts/compare_comsol/compare_comsol_DFN.py
+++ b/examples/scripts/compare_comsol/compare_comsol_DFN.py
@@ -2,8 +2,7 @@ import pybamm
 import numpy as np
 import os
 import pickle
-import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider
+import scipy.interpolate as interp
 
 # change working directory to the root of pybamm
 os.chdir(pybamm.__path__[0] + "/..")
@@ -26,234 +25,93 @@ comsol_variables = pickle.load(
 "Create and solve pybamm model"
 
 # load model and geometry
-model = pybamm.lithium_ion.DFN()
-geometry = model.default_geometry
+pybamm.set_logging_level("INFO")
+pybamm_model = pybamm.lithium_ion.DFN()
+geometry = pybamm_model.default_geometry
 
 # load parameters and process model and geometry
-param = model.default_parameter_values
+param = pybamm_model.default_parameter_values
 param["Electrode depth [m]"] = 1
 param["Electrode height [m]"] = 1
 param["Typical current [A]"] = 24 * C_rates[C_rate]
-param.process_model(model)
+param.process_model(pybamm_model)
 param.process_geometry(geometry)
 
 # create mesh
 var = pybamm.standard_spatial_vars
 var_pts = {var.x_n: 31, var.x_s: 11, var.x_p: 31, var.r_n: 11, var.r_p: 11}
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+mesh = pybamm.Mesh(geometry, pybamm_model.default_submesh_types, var_pts)
 
 # discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
+disc = pybamm.Discretisation(mesh, pybamm_model.default_spatial_methods)
+disc.process_model(pybamm_model)
 
 # discharge timescale
-tau = param.process_symbol(
-    pybamm.standard_parameters_lithium_ion.tau_discharge
-).evaluate(0, 0)
+tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge)
 
 # solve model at comsol times
-time = comsol_variables["time"] / tau
-solution = model.default_solver.solve(model, time)
-
-"-----------------------------------------------------------------------------"
-"Get variables for comparison"
-
-# discharge capacity
-discharge_capacity = pybamm.ProcessedVariable(
-    model.variables["Discharge capacity [A.h]"], solution.t, solution.y, mesh=mesh
-)
-discharge_capacity_sol = discharge_capacity(solution.t)
-comsol_discharge_capacity = (
-    comsol_variables["time"] * param["Typical current [A]"] / 3600
-)
-
-# spatial points
-l_n = param.process_symbol(pybamm.geometric_parameters.l_n).evaluate(0, 0)
-l_s = param.process_symbol(pybamm.geometric_parameters.l_s).evaluate(0, 0)
-l_p = param.process_symbol(pybamm.geometric_parameters.l_p).evaluate(0, 0)
-x_n = np.linspace(0, l_n, 40)
-x_s = np.linspace(l_n, l_n + l_s, 20)
-x_p = np.linspace(l_n + l_s, 1, 40)
-x = np.linspace(0, 1, 100)
-
-# output variables
-output_variables = {
-    "c_n_surf": model.variables["Negative particle surface concentration [mol.m-3]"],
-    "c_e": model.variables["Electrolyte concentration [mol.m-3]"],
-    "c_p_surf": model.variables["Positive particle surface concentration [mol.m-3]"],
-    "phi_n": model.variables["Negative electrode potential [V]"],
-    "phi_e": model.variables["Electrolyte potential [V]"],
-    "phi_p": model.variables["Positive electrode potential [V]"],
-    "voltage": model.variables["Terminal voltage [V]"],
+time = comsol_variables["time"] / tau.evaluate(0)
+solution = pybamm_model.default_solver.solve(pybamm_model, time)
+
+
+# Make Comsol 'model' for comparison
+whole_cell = ["negative electrode", "separator", "positive electrode"]
+comsol_t = comsol_variables["time"]
+L_x = param.process_symbol(pybamm.standard_parameters_lithium_ion.L_x).evaluate()
+
+
+def get_interp_fun(variable, domain):
+    """
+    Create a :class:`pybamm.Function` object using the variable, to allow plotting with
+    :class:`'pybamm.QuickPlot'` (interpolate in space to match edges, and then create
+    function to interpolate in time)
+    """
+    if domain == ["negative electrode"]:
+        comsol_x = comsol_variables["x_n"]
+    elif domain == ["positive electrode"]:
+        comsol_x = comsol_variables["x_p"]
+    elif domain == whole_cell:
+        comsol_x = comsol_variables["x"]
+    # Make sure to use dimensional space
+    pybamm_x = mesh.combine_submeshes(*domain)[0].nodes * L_x
+    variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x)
+
+    def myinterp(t):
+        return interp.interp1d(comsol_t, variable)(t)[:, np.newaxis]
+
+    # Make sure to use dimensional time
+    fun = pybamm.Function(myinterp, pybamm.t * tau)
+    fun.domain = domain
+    return fun
+
+
+comsol_c_n_surf = get_interp_fun(comsol_variables["c_n_surf"], ["negative electrode"])
+comsol_c_e = get_interp_fun(comsol_variables["c_e"], whole_cell)
+comsol_c_p_surf = get_interp_fun(comsol_variables["c_p_surf"], ["positive electrode"])
+comsol_phi_n = get_interp_fun(comsol_variables["phi_n"], ["negative electrode"])
+comsol_phi_e = get_interp_fun(comsol_variables["phi_e"], whole_cell)
+comsol_phi_p = get_interp_fun(comsol_variables["phi_p"], ["positive electrode"])
+comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"])
+
+# Create comsol model with dictionary of Matrix variables
+comsol_model = pybamm.BaseModel()
+comsol_model.variables = {
+    "Negative particle surface concentration [mol.m-3]": comsol_c_n_surf,
+    "Electrolyte concentration [mol.m-3]": comsol_c_e,
+    "Positive particle surface concentration [mol.m-3]": comsol_c_p_surf,
+    "Current [A]": pybamm_model.variables["Current [A]"],
+    "Negative electrode potential [V]": comsol_phi_n,
+    "Electrolyte potential [V]": comsol_phi_e,
+    "Positive electrode potential [V]": comsol_phi_p,
+    "Terminal voltage [V]": pybamm.Function(comsol_voltage, pybamm.t * tau),
 }
 
-processed_variables = pybamm.post_process_variables(
-    output_variables, solution.t, solution.y, mesh=mesh
-)
-
-"-----------------------------------------------------------------------------"
-"Make plots"
-
-fig, ax = plt.subplots(figsize=(15, 8))
-plt.tight_layout()
-plt.subplots_adjust(left=-0.1)
-
-
-# negative particle surface concentration
-c_n_surf_min = 0.9 * np.min(comsol_variables["c_n_surf"])
-c_n_surf_max = 1.1 * np.max(comsol_variables["c_n_surf"])
-plt.subplot(241)
-c_n_surf_plot, = plt.plot(x_n, processed_variables["c_n_surf"](time[0], x_n), "b-")
-comsol_c_n_surf_plot, = plt.plot(
-    comsol_variables["x_n"] / comsol_variables["x"][-1],
-    comsol_variables["c_n_surf"][:, 0],
-    "r:",
-)
-plt.axis([0, l_n, c_n_surf_min, c_n_surf_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"Surface $c_n$ (mol/m$^3$)")
-
-# electrolyte concentration
-c_e_min = 0.9 * np.min(comsol_variables["c_e"])
-c_e_max = 1.1 * np.max(comsol_variables["c_e"])
-plt.subplot(242)
-c_e_plot, = plt.plot(x, processed_variables["c_e"](time[0], x), "b-")
-comsol_c_e_plot, = plt.plot(
-    comsol_variables["x"] / comsol_variables["x"][-1],
-    comsol_variables["c_e"][:, 0],
-    "r:",
-)
-plt.axis([0, 1, c_e_min, c_e_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$c_e$ (mol/m$^3$)")
-
-# negative particle surface concentration
-c_p_surf_min = 0.9 * np.min(comsol_variables["c_p_surf"])
-c_p_surf_max = 1.1 * np.max(comsol_variables["c_p_surf"])
-plt.subplot(243)
-c_p_surf_plot, = plt.plot(x_p, processed_variables["c_n_surf"](time[0], x_p), "b-")
-comsol_c_p_surf_plot, = plt.plot(
-    comsol_variables["x_p"] / comsol_variables["x"][-1],
-    comsol_variables["c_p_surf"][:, 0],
-    "r:",
-)
-plt.axis([l_n + l_s, 1, c_p_surf_min, c_p_surf_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"Surface $c_p$ (mol/m$^3$)")
-
-# current
-I_min = 0.9 * param["Typical current [A]"]
-I_max = 1.1 * param["Typical current [A]"]
-plt.subplot(244)
-time_tracer, = plt.plot(
-    [comsol_variables["time"][0], comsol_variables["time"][0]], [I_min, I_max], "k--"
-)
-plt.plot(
-    [comsol_variables["time"][0], comsol_variables["time"][-1]],
-    [param["Typical current [A]"], param["Typical current [A]"]],
-    "b-",
-)
-plt.axis(
-    [comsol_variables["time"][0], comsol_variables["time"][-1] * 1.1, I_min, I_max]
-)
-plt.xlabel("Time (s)")
-plt.ylabel("Applied current (A)")
-
-# negative electrode potential
-phi_n_min = 1.1 * np.min(comsol_variables["phi_n"])
-phi_n_max = 0.9 * np.max(comsol_variables["phi_n"])
-plt.subplot(245)
-phi_n_plot, = plt.plot(x_n, processed_variables["phi_n"](time[0], x_n), "b-")
-comsol_phi_n_plot, = plt.plot(
-    comsol_variables["x_n"] / comsol_variables["x"][-1],
-    comsol_variables["phi_n"][:, 0],
-    "r:",
-)
-plt.axis([0, l_n, phi_n_min, phi_n_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$\phi_n$ (V)")
-
-# electrolyte potential
-phi_e_min = 1.1 * np.min(comsol_variables["phi_e"])
-phi_e_max = 0.9 * np.max(comsol_variables["phi_e"])
-plt.subplot(246)
-phi_e_plot, = plt.plot(x, processed_variables["phi_e"](time[0], x), "b-")
-comsol_phi_e_plot, = plt.plot(
-    comsol_variables["x"] / comsol_variables["x"][-1],
-    comsol_variables["phi_e"][:, 0],
-    "r:",
-)
-plt.axis([0, 1, phi_e_min, phi_e_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$\phi_e$ (V)")
-
-# positive electrode potential
-phi_p_min = 0.9 * np.min(comsol_variables["phi_p"])
-phi_p_max = 1.1 * np.max(comsol_variables["phi_p"])
-plt.subplot(247)
-phi_p_plot, = plt.plot(x_p, processed_variables["phi_p"](time[0], x_p), "b-")
-comsol_phi_p_plot, = plt.plot(
-    comsol_variables["x_p"] / comsol_variables["x"][-1],
-    comsol_variables["phi_p"][:, 0],
-    "r:",
+# plot
+plot = pybamm.QuickPlot(
+    [pybamm_model, comsol_model],
+    mesh,
+    [solution, solution],
+    output_variables=comsol_model.variables.keys(),
+    labels=["PyBaMM", "Comsol"],
 )
-plt.axis([l_n + l_s, 1, phi_p_min, phi_p_max])
-plt.xlabel(r"$x$")
-plt.ylabel(r"$\phi_p$ (V)")
-
-# discharge curve
-v_min = 3.2
-v_max = 3.9
-plt.subplot(248)
-discharge_capacity_tracer, = plt.plot(
-    [comsol_discharge_capacity[0], comsol_discharge_capacity[0]], [v_min, v_max], "k--"
-)
-plt.plot(
-    discharge_capacity_sol,
-    processed_variables["voltage"](solution.t),
-    "b-",
-    label="PyBaMM",
-)
-plt.plot(comsol_discharge_capacity, comsol_variables["voltage"], "r:", label="Comsol")
-plt.axis(
-    [comsol_discharge_capacity[0], comsol_discharge_capacity[-1] * 1.1, v_min, v_max]
-)
-plt.xlabel("Discharge Capacity (Ah)")
-plt.ylabel("Voltage (V)")
-plt.legend(loc="best")
-
-axcolor = "lightgoldenrodyellow"
-axfreq = plt.axes([0.315, 0.02, 0.37, 0.03], facecolor=axcolor)
-sfreq = Slider(axfreq, "Time", 0, comsol_variables["time"][-1], valinit=0)
-
-
-def update_plot(t):
-    # find t index
-    ind = (np.abs(comsol_variables["time"] - t)).argmin()
-    # update time
-    time_tracer.set_xdata(comsol_variables["time"][ind])
-    discharge_capacity_tracer.set_xdata(comsol_discharge_capacity[ind])
-    # update negative particle surface concentration
-    c_n_surf_plot.set_ydata(processed_variables["c_n_surf"](time[ind], x_n))
-    comsol_c_n_surf_plot.set_ydata(comsol_variables["c_n_surf"][:, ind])
-    # update electrolyte concentration
-    c_e_plot.set_ydata(processed_variables["c_e"](time[ind], x))
-    comsol_c_e_plot.set_ydata(comsol_variables["c_e"][:, ind])
-    # update negative particle surface concentration
-    c_p_surf_plot.set_ydata(processed_variables["c_p_surf"](time[ind], x_p))
-    comsol_c_p_surf_plot.set_ydata(comsol_variables["c_p_surf"][:, ind])
-    # update negative electrode potential
-    phi_n_plot.set_ydata(processed_variables["phi_n"](time[ind], x_n))
-    comsol_phi_n_plot.set_ydata(comsol_variables["phi_n"][:, ind])
-    # update electrolyte potential
-    phi_e_plot.set_ydata(processed_variables["phi_e"](time[ind], x))
-    comsol_phi_e_plot.set_ydata(comsol_variables["phi_e"][:, ind])
-    # update positive electrode potential
-    phi_p_plot.set_ydata(processed_variables["phi_p"](time[ind], x_p))
-    comsol_phi_p_plot.set_ydata(comsol_variables["phi_p"][:, ind])
-    fig.canvas.draw_idle()
-
-
-sfreq.on_changed(update_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.dynamic_plot()
diff --git a/examples/scripts/compare_lead_acid.py b/examples/scripts/compare_lead_acid.py
index cdf22b3db..4ecccb998 100644
--- a/examples/scripts/compare_lead_acid.py
+++ b/examples/scripts/compare_lead_acid.py
@@ -47,7 +47,10 @@ for i, model in enumerate(models):
 
 # plot
 output_variables = [
-    "Interfacial current density [A.m-2]",
+    [
+        "Average negative electrode interfacial current density [A.m-2]",
+        "Average positive electrode interfacial current density [A.m-2]",
+    ],
     "Electrolyte concentration [mol.m-3]",
     "Porosity",
     "Electrolyte current density [A.m-2]",
diff --git a/pandas/_libs/parsers.pyx b/pandas/_libs/parsers.pyx
deleted file mode 100644
index e69de29bb..000000000
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 8701785c6..efc634fa6 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -25,7 +25,10 @@ class Function(pybamm.Symbol):
 
     def __init__(self, function, *children):
 
-        name = "function ({})".format(function.__name__)
+        try:
+            name = "function ({})".format(function.__name__)
+        except AttributeError:
+            name = "function ({})".format(function.__class__)
         children_list = list(children)
         domain = self.get_children_domains(children_list)
 
diff --git a/pybamm/quick_plot.py b/pybamm/quick_plot.py
index 79a23a55d..5465d6be7 100644
--- a/pybamm/quick_plot.py
+++ b/pybamm/quick_plot.py
@@ -3,6 +3,7 @@
 #
 import numpy as np
 import pybamm
+from collections import defaultdict
 
 
 def ax_min(data):
@@ -23,6 +24,18 @@ def ax_max(data):
         return 1.1 * data_max
 
 
+def split_long_string(title, max_words=4):
+    "Get title in a nice format"
+    words = title.split()
+    # Don't split if fits on one line, don't split just for units
+    if len(words) <= max_words or words[max_words].startswith("["):
+        return title
+    else:
+        first_line = (" ").join(words[:max_words])
+        second_line = (" ").join(words[max_words:])
+        return first_line + "\n" + second_line
+
+
 class QuickPlot(object):
     """
     Generates a quick plot of a subset of key outputs of the model so that the model
@@ -39,11 +52,14 @@ class QuickPlot(object):
     solutions: (iter of) :class:`pybamm.Solver`
         The numerical solution(s) for the model(s) which contained the solution to the
         model(s).
-    output_variables : list of str
+    output_variables : list of str, optional
         List of variables to plot
+    labels : list of str, optional
+        Labels for the different models. Defaults to model names
     """
 
-    def __init__(self, models, mesh, solutions, output_variables=None):
+    def __init__(self, models, mesh, solutions, output_variables=None, labels=None):
+        # Pre-process models and solutions
         if isinstance(models, pybamm.BaseModel):
             models = [models]
         elif not isinstance(models, list):
@@ -57,14 +73,17 @@ class QuickPlot(object):
         else:
             raise ValueError("must provide the same number of models and solutions")
 
+        # Set labels
+        self.labels = labels or [model.name for model in models]
+
         # Scales (default to 1 if information not in model)
-        vars = models[0].variables
+        variables = models[0].variables
         self.x_scale = 1
         self.time_scale = 1
-        if "x [m]" and "x" in vars:
-            self.x_scale = (vars["x [m]"] / vars["x"]).evaluate()[-1]
-        if "Time [m]" and "Time" in vars:
-            self.time_scale = (vars["Time [h]"] / vars["Time"]).evaluate(t=1)
+        if "x [m]" and "x" in variables:
+            self.x_scale = (variables["x [m]"] / variables["x"]).evaluate()[-1]
+        if "Time [m]" and "Time" in variables:
+            self.time_scale = (variables["Time [h]"] / variables["Time"]).evaluate(t=1)
 
         # Time parameters
         self.ts = [solution.t for solution in solutions]
@@ -107,55 +126,83 @@ class QuickPlot(object):
 
         # Calculate subplot positions based on number of variables supplied
         self.subplot_positions = {}
-        n = int(len(output_variables) // np.sqrt(len(output_variables)))
-        m = np.ceil(len(output_variables) / n)
+        self.n_rows = int(len(output_variables) // np.sqrt(len(output_variables)))
+        self.n_cols = int(np.ceil(len(output_variables) / self.n_rows))
 
         # Process output variables into a form that can be plotted
-        for k, var in enumerate(output_variables):
-            self.variables[var] = [
-                pybamm.ProcessedVariable(
-                    models[i].variables[var], solutions[i].t, solutions[i].y, mesh
-                )
-                for i in range(len(models))
-            ]
-            if self.variables[var][0].dimensions == 2:
-                domain = models[0].variables[var].domain
-                self.x_values[var] = mesh.combine_submeshes(*domain)[0].edges
-            self.subplot_positions[var] = (n, m, k + 1)
+        for k, variable_list in enumerate(output_variables):
+            # Make sure we always have a list of lists of variables
+            if isinstance(variable_list, str):
+                variable_list = [variable_list]
 
-        # Set labels
-        self.labels = [model.name for model in models]
+            # Prepare list of variables
+            key = tuple(variable_list)
+            self.variables[key] = [None] * len(models)
+
+            # process each variable in variable_list for each model
+            for i, model in enumerate(models):
+                # self.variables is a dictionary of lists of lists
+                self.variables[key][i] = [
+                    pybamm.ProcessedVariable(
+                        model.variables[var], solutions[i].t, solutions[i].y, mesh
+                    )
+                    for var in variable_list
+                ]
+
+            # Make sure variables have the same dimensions and domain
+            domain = self.variables[key][0][0].domain
+            for variable in self.variables[key][0]:
+                if variable.domain != domain:
+                    raise ValueError("mismatching variable domains")
+
+            # Set the x variable for any two-dimensional variables
+            if self.variables[key][0][0].dimensions == 2:
+                self.x_values[key] = mesh.combine_submeshes(*domain)[0].edges
+
+            # Don't allow 3D variables
+            elif any(var.dimensions == 3 for var in self.variables[key][0]):
+                raise NotImplementedError("cannot plot 3D variables")
 
-        # Don't allow 3D variables
-        if self.variables[var][0].dimensions == 3:
-            raise NotImplementedError("cannot plot 3D variables")
+            # Define subplot position
+            self.subplot_positions[key] = (self.n_rows, self.n_cols, k + 1)
 
     def reset_axis(self):
         """
         Reset the axis limits to the default values.
+        These are calculated to fit around the minimum and maximum values of all the
+        variables in each subplot
         """
         self.axis = {}
-        for name, variable in self.variables.items():
-            if variable[0].dimensions == 1:
-                y_min = np.min([ax_min(v(self.ts[i])) for i, v in enumerate(variable)])
-                y_max = np.max([ax_max(v(self.ts[i])) for i, v in enumerate(variable)])
-                if y_min == y_max:
-                    y_min -= 1
-                    y_max += 1
-                self.axis[name] = [self.min_t, self.max_t, y_min, y_max]
-            elif variable[0].dimensions == 2:
-                x = self.x_values[name]
+        for key, variable_lists in self.variables.items():
+            if variable_lists[0][0].dimensions == 1:
+                x = None
+                x_min = self.min_t
+                x_max = self.max_t
+            elif variable_lists[0][0].dimensions == 2:
+                x = self.x_values[key]
                 x_scaled = x * self.x_scale
-                y_min = np.min(
-                    [ax_min(v(self.ts[i], x)) for i, v in enumerate(variable)]
-                )
-                y_max = np.max(
-                    [ax_max(v(self.ts[i], x)) for i, v in enumerate(variable)]
-                )
-                if y_min == y_max:
-                    y_min -= 1
-                    y_max += 1
-                self.axis[name] = [x_scaled[0], x_scaled[-1], y_min, y_max]
+                x_min = x_scaled[0]
+                x_max = x_scaled[-1]
+
+            # Get min and max y values
+            y_min = np.min(
+                [
+                    ax_min(var(self.ts[i], x))
+                    for i, variable_list in enumerate(variable_lists)
+                    for var in variable_list
+                ]
+            )
+            y_max = np.max(
+                [
+                    ax_max(var(self.ts[i], x))
+                    for i, variable_list in enumerate(variable_lists)
+                    for var in variable_list
+                ]
+            )
+            if y_min == y_max:
+                y_min -= 1
+                y_max += 1
+            self.axis[key] = [x_min, x_max, y_min, y_max]
 
     def plot(self, t):
         """Produces a quick plot with the internal states at time t.
@@ -169,51 +216,81 @@ class QuickPlot(object):
         import matplotlib.pyplot as plt
 
         t /= self.time_scale
-        self.fig, self.ax = plt.subplots(figsize=(15, 8))
+        self.fig, self.ax = plt.subplots(self.n_rows, self.n_cols, figsize=(15, 8))
         plt.tight_layout()
         plt.subplots_adjust(left=-0.1)
         self.plots = {}
         self.time_lines = {}
 
-        for k, name in enumerate(self.variables.keys()):
-            variable = self.variables[name]
-            plt.subplot(*self.subplot_positions[name])
-            plt.ylabel(name, fontsize=14)
-            plt.axis(self.axis[name])
-            self.plots[name] = [None] * self.num_models
-            # Set labels
+        colors = ["r", "b", "k", "g"]
+        linestyles = ["-", ":", "--", "-."]
+        fontsize = 42 // self.n_cols
+
+        for k, (key, variable_lists) in enumerate(self.variables.items()):
+            if len(self.variables) == 1:
+                ax = self.ax
+            else:
+                ax = self.ax.flat[k]
+            ax.set_xlim(self.axis[key][:2])
+            ax.set_ylim(self.axis[key][2:])
+            ax.xaxis.set_major_locator(plt.MaxNLocator(3))
+            self.plots[key] = defaultdict(dict)
+            # Set labels for the first subplot only (avoid repetition)
             if k == 0:
                 labels = self.labels
             else:
                 labels = [None] * len(self.labels)
-            if variable[0].dimensions == 2:
+            if variable_lists[0][0].dimensions == 2:
                 # 2D plot: plot as a function of x at time t
-                plt.xlabel("Position [m]", fontsize=14)
-                x_value = self.x_values[name]
-                for i in range(self.num_models):
-                    self.plots[name][i], = plt.plot(
-                        x_value * self.x_scale,
-                        variable[i](t, x_value),
-                        lw=2,
-                        label=labels[i],
-                    )
+                ax.set_xlabel("Position [m]", fontsize=fontsize)
+                x_value = self.x_values[key]
+                for i, variable_list in enumerate(variable_lists):
+                    for j, variable in enumerate(variable_list):
+                        if j == 0:
+                            label = labels[i]
+                        else:
+                            label = None
+                        self.plots[key][i][j], = ax.plot(
+                            x_value * self.x_scale,
+                            variable(t, x_value),
+                            lw=2,
+                            color=colors[i],
+                            linestyle=linestyles[j],
+                            label=label,
+                        )
             else:
                 # 1D plot: plot as a function of time, indicating time t with a line
-                plt.xlabel("Time [h]", fontsize=14)
-                for i in range(self.num_models):
-                    full_t = self.ts[i]
-                    self.plots[name][i], = plt.plot(
-                        full_t * self.time_scale,
-                        variable[i](full_t),
-                        lw=2,
-                        label=labels[i],
-                    )
-                    y_min, y_max = self.axis[name][2:]
-                    self.time_lines[name], = plt.plot(
-                        [t * self.time_scale, t * self.time_scale],
-                        [y_min, y_max],
-                        "k--",
-                    )
+                ax.set_xlabel("Time [h]", fontsize=fontsize)
+                for i, variable_list in enumerate(variable_lists):
+                    for j, variable in enumerate(variable_list):
+                        full_t = self.ts[i]
+                        if j == 0:
+                            label = labels[i]
+                        else:
+                            label = None
+                        self.plots[key][i][j], = ax.plot(
+                            full_t * self.time_scale,
+                            variable(full_t),
+                            lw=2,
+                            color=colors[i],
+                            linestyle=linestyles[j],
+                            label=label,
+                        )
+                y_min, y_max = self.axis[key][2:]
+                self.time_lines[key], = ax.plot(
+                    [t * self.time_scale, t * self.time_scale], [y_min, y_max], "k--"
+                )
+            # Set either y label or legend entries
+            if len(key) == 1:
+                title = split_long_string(key[0])
+                ax.set_title(title, fontsize=fontsize)
+            else:
+                ax.legend(
+                    [split_long_string(s, 6) for s in key],
+                    bbox_to_anchor=(0.5, 1.2),
+                    fontsize=8,
+                    loc="upper center",
+                )
         self.fig.legend(loc="lower right")
 
     def dynamic_plot(self, testing=False):
@@ -246,12 +323,13 @@ class QuickPlot(object):
         """
         t = self.sfreq.val
         t_dimensionless = t / self.time_scale
-        for var, plot in self.plots.items():
-            if self.variables[var][0].dimensions == 2:
-                x = self.x_values[var]
-                for i in range(self.num_models):
-                    plot[i].set_ydata(self.variables[var][i](t_dimensionless, x))
+        for key, plot in self.plots.items():
+            if self.variables[key][0][0].dimensions == 2:
+                x = self.x_values[key]
+                for i, variable_lists in enumerate(self.variables[key]):
+                    for j, variable in enumerate(variable_lists):
+                        plot[i][j].set_ydata(variable(t_dimensionless, x))
             else:
-                self.time_lines[var].set_xdata([t])
+                self.time_lines[key].set_xdata([t])
 
         self.fig.canvas.draw_idle()

Quick Plot functionality
Add more functionality to `QuickPlot`:
- [x] User can define labels (otherwise defaults to model names)
- [x] Can provide list of list of variables, which plots several variables on the same graph
- [ ] Adjust x ticks to fit for lithium ion plots
- [x] Automatically adjust label sizes to avoid overlap
diff --git a/tests/integration/test_quick_plot.py b/tests/integration/test_quick_plot.py
index 05b49d33e..d573bf674 100644
--- a/tests/integration/test_quick_plot.py
+++ b/tests/integration/test_quick_plot.py
@@ -30,12 +30,12 @@ class TestQuickPlot(unittest.TestCase):
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        quick_plot.axis.update({"Electrolyte concentration": new_axis})
-        self.assertEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        quick_plot.axis.update({("Electrolyte concentration",): new_axis})
+        self.assertEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        self.assertNotEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
@@ -61,12 +61,12 @@ class TestQuickPlot(unittest.TestCase):
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        quick_plot.axis.update({"Electrolyte concentration": new_axis})
-        self.assertEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        quick_plot.axis.update({("Electrolyte concentration",): new_axis})
+        self.assertEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis["Electrolyte concentration"], new_axis)
+        self.assertNotEqual(quick_plot.axis[("Electrolyte concentration",)], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
diff --git a/tests/unit/test_expression_tree/test_functions.py b/tests/unit/test_expression_tree/test_functions.py
index 70b1e8a04..0ac27a800 100644
--- a/tests/unit/test_expression_tree/test_functions.py
+++ b/tests/unit/test_expression_tree/test_functions.py
@@ -6,6 +6,7 @@ import pybamm
 import unittest
 import numpy as np
 import autograd.numpy as auto_np
+from scipy.interpolate import interp1d
 
 
 def test_function(arg):
@@ -93,6 +94,15 @@ class TestFunction(unittest.TestCase):
         with self.assertRaises(pybamm.DomainError):
             pybamm.Function(test_multi_var_function, a, b)
 
+    def test_function_unnamed(self):
+        t = np.linspace(0, 1)
+        entries = 2 * t
+        interpfun = interp1d(t, entries)
+        fun = pybamm.Function(interpfun, pybamm.t)
+        self.assertEqual(
+            fun.name, "function (<class 'scipy.interpolate.interpolate.interp1d'>)"
+        )
+
 
 class TestSpecificFunctions(unittest.TestCase):
     def test_cos(self):
diff --git a/tests/unit/test_quick_plot.py b/tests/unit/test_quick_plot.py
index 526a04af4..3913c50f2 100644
--- a/tests/unit/test_quick_plot.py
+++ b/tests/unit/test_quick_plot.py
@@ -21,12 +21,12 @@ class TestQuickPlot(unittest.TestCase):
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        quick_plot.axis.update({"a": new_axis})
-        self.assertEqual(quick_plot.axis["a"], new_axis)
+        quick_plot.axis.update({("a",): new_axis})
+        self.assertEqual(quick_plot.axis[("a",)], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis["a"], new_axis)
+        self.assertNotEqual(quick_plot.axis[("a",)], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
@@ -34,27 +34,66 @@ class TestQuickPlot(unittest.TestCase):
         quick_plot.update(0.01)
 
         # Test with different output variables
+        quick_plot = pybamm.QuickPlot(model, mesh, solution, ["b broadcasted"])
+        self.assertEqual(len(quick_plot.axis), 1)
+        quick_plot.plot(0)
+
         quick_plot = pybamm.QuickPlot(
-            model, mesh, solution, ["b broadcasted", "c broadcasted"]
+            model,
+            mesh,
+            solution,
+            [["a", "a"], ["b broadcasted", "b broadcasted"], "c broadcasted"],
         )
-        self.assertEqual(len(quick_plot.axis), 2)
+        self.assertEqual(len(quick_plot.axis), 3)
         quick_plot.plot(0)
 
         # update the axis
         new_axis = [0, 0.5, 0, 1]
-        var = "c broadcasted"
-        quick_plot.axis.update({var: new_axis})
-        self.assertEqual(quick_plot.axis[var], new_axis)
+        var_key = ("c broadcasted",)
+        quick_plot.axis.update({var_key: new_axis})
+        self.assertEqual(quick_plot.axis[var_key], new_axis)
 
         # and now reset them
         quick_plot.reset_axis()
-        self.assertNotEqual(quick_plot.axis[var], new_axis)
+        self.assertNotEqual(quick_plot.axis[var_key], new_axis)
 
         # check dynamic plot loads
         quick_plot.dynamic_plot(testing=True)
 
         quick_plot.update(0.01)
 
+        # Test longer name
+        model.variables["Variable with a very long name"] = model.variables["a"]
+        quick_plot = pybamm.QuickPlot(model, mesh, solution)
+        quick_plot.plot(0)
+
+        # Test errors
+        with self.assertRaisesRegex(ValueError, "mismatching variable domains"):
+            pybamm.QuickPlot(model, mesh, solution, [["a", "b broadcasted"]])
+        model.variables["3D variable"] = disc.process_symbol(
+            pybamm.Broadcast(1, ["negative particle"])
+        )
+        with self.assertRaisesRegex(NotImplementedError, "cannot plot 3D variables"):
+            pybamm.QuickPlot(model, mesh, solution, ["3D variable"])
+
+    def test_loqs_spm_base(self):
+        t_eval = np.linspace(0, 0.01, 2)
+
+        # SPM
+        for model in [pybamm.lithium_ion.SPM(), pybamm.lead_acid.LOQS()]:
+            geometry = model.default_geometry
+            param = model.default_parameter_values
+            param.process_model(model)
+            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)
+            solver = model.default_solver
+            solution = solver.solve(model, t_eval)
+            pybamm.QuickPlot(model, mesh, solution)
+
     def test_failure(self):
         with self.assertRaisesRegex(TypeError, "'models' must be"):
             pybamm.QuickPlot(1, None, None)

pip install -e .[dev,docs]
parse_log_pytest
pytest
pytest
apt-get update apt-get install -y gfortran gcc libopenblas-dev liblapack-dev
3.6
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs @ file:///opt/conda/conda-bld/attrs_1642510447205/work
autograd==1.6.2
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
flake8==5.0.4
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
more-itertools @ file:///tmp/build/80754af9/more-itertools_1637733554872/work
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging @ file:///tmp/build/80754af9/packaging_1637314298585/work
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy @ file:///tmp/build/80754af9/pluggy_1615976315926/work
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py @ file:///opt/conda/conda-bld/py_1644396412707/work
-e git+https://github.com/pybamm-team/PyBaMM.git@4f87f454d8266f46553f304a4f9895bd1720d438#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing @ file:///tmp/build/80754af9/pyparsing_1635766073266/work
pyrsistent==0.18.0
pytest==6.2.4
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
toml @ file:///tmp/build/80754af9/toml_1616166611790/work
tornado==6.1
traitlets==4.3.3
typing_extensions @ file:///opt/conda/conda-bld/typing_extensions_1647553014482/work
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp @ file:///tmp/build/80754af9/zipp_1633618647012/work

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - attrs=21.4.0=pyhd3eb1b0_0
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - importlib_metadata=4.8.1=hd3eb1b0_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - more-itertools=8.12.0=pyhd3eb1b0_0
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - packaging=21.3=pyhd3eb1b0_0
  - pip=21.2.2=py36h06a4308_0
  - pluggy=0.13.1=py36h06a4308_0
  - py=1.11.0=pyhd3eb1b0_0
  - pyparsing=3.0.4=pyhd3eb1b0_0
  - pytest=6.2.4=py36h06a4308_2
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - toml=0.10.2=pyhd3eb1b0_0
  - typing_extensions=4.1.1=pyh06a4308_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zipp=3.6.0=pyhd3eb1b0_0
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - autograd==1.6.2
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - flake8==5.0.4
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyrsistent==0.18.0
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tornado==6.1
      - traitlets==4.3.3
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_expression_tree/test_functions.py::TestFunction::test_function_unnamed tests/unit/test_quick_plot.py::TestQuickPlot::test_simple_ode_model
tests/unit/test_expression_tree/test_functions.py::test_const_function tests/unit/test_expression_tree/test_functions.py::TestFunction::test_constant_functions tests/unit/test_expression_tree/test_functions.py::TestFunction::test_exceptions tests/unit/test_expression_tree/test_functions.py::TestFunction::test_function_of_multiple_variables tests/unit/test_expression_tree/test_functions.py::TestFunction::test_function_of_one_variable tests/unit/test_expression_tree/test_functions.py::TestFunction::test_with_autograd tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_cos tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_cosh tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_exp tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_log tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_max tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_min tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_sin tests/unit/test_expression_tree/test_functions.py::TestSpecificFunctions::test_sinh tests/unit/test_quick_plot.py::TestQuickPlot::test_failure tests/unit/test_quick_plot.py::TestQuickPlot::test_loqs_spm_base
examples/scripts/compare_lead_acid.py pybamm/expression_tree/functions.py pandas/_libs/parsers.pyx pybamm/quick_plot.py .gitignore examples/scripts/compare_comsol/compare_comsol_DFN.py
examples/scripts/compare_lead_acid.py pybamm/expression_tree/functions.py pandas/_libs/parsers.pyx pybamm/quick_plot.py .gitignore examples/scripts/compare_comsol/compare_comsol_DFN.py