diff --git a/CHANGELOG.md b/CHANGELOG.md
index f7e08058c..6bce9313f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+- [#450](https://github.com/pybop-team/PyBOP/pull/450) - Adds support for IDAKLU with output variables, and corresponding examples, tests.
 - [#364](https://github.com/pybop-team/PyBOP/pull/364) - Adds the MultiFittingProblem class and the multi_fitting example script.
 - [#444](https://github.com/pybop-team/PyBOP/issues/444) - Merge `BaseModel` `build()` and `rebuild()` functionality.
 - [#435](https://github.com/pybop-team/PyBOP/pull/435) - Adds SLF001 linting for private members.
diff --git a/examples/notebooks/equivalent_circuit_identification.ipynb b/examples/notebooks/equivalent_circuit_identification.ipynb
index 6aa0d8907..1e7b10b92 100644
--- a/examples/notebooks/equivalent_circuit_identification.ipynb
+++ b/examples/notebooks/equivalent_circuit_identification.ipynb
@@ -11,7 +11,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/notebooks/multi_model_identification.ipynb b/examples/notebooks/multi_model_identification.ipynb
index 1b26b2862..e9b3e9b53 100644
--- a/examples/notebooks/multi_model_identification.ipynb
+++ b/examples/notebooks/multi_model_identification.ipynb
@@ -12,7 +12,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/notebooks/multi_optimiser_identification.ipynb b/examples/notebooks/multi_optimiser_identification.ipynb
index 2131c2b38..7067f6bf0 100644
--- a/examples/notebooks/multi_optimiser_identification.ipynb
+++ b/examples/notebooks/multi_optimiser_identification.ipynb
@@ -12,7 +12,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/notebooks/optimiser_calibration.ipynb b/examples/notebooks/optimiser_calibration.ipynb
index 2b0fd3348..d1c67b08e 100644
--- a/examples/notebooks/optimiser_calibration.ipynb
+++ b/examples/notebooks/optimiser_calibration.ipynb
@@ -12,7 +12,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/notebooks/pouch_cell_identification.ipynb b/examples/notebooks/pouch_cell_identification.ipynb
index 146ee64bd..2f52ab825 100644
--- a/examples/notebooks/pouch_cell_identification.ipynb
+++ b/examples/notebooks/pouch_cell_identification.ipynb
@@ -12,7 +12,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/notebooks/solver_selection.ipynb b/examples/notebooks/solver_selection.ipynb
new file mode 100644
index 000000000..81db4b6c9
--- /dev/null
+++ b/examples/notebooks/solver_selection.ipynb
@@ -0,0 +1,261 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "id": "expmkveO04pw"
+   },
+   "source": [
+    "## Investigate Different PyBaMM Solvers\n",
+    "\n",
+    "In this notebook, we discuss the process of changing PyBaMM solvers and the corresponding performance trade-offs with each. For further reading on different solvers, see the PyBaMM solver documentation:\n",
+    "\n",
+    "[[1]: PyBaMM Solvers](https://docs.pybamm.org/en/stable/source/api/solvers/index.html#)\n",
+    "\n",
+    "### Setting up the Environment\n",
+    "\n",
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "colab": {
+     "base_uri": "https://localhost:8080/"
+    },
+    "id": "X87NUGPW04py",
+    "outputId": "0d785b07-7cff-4aeb-e60a-4ff5a669afbf"
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Note: you may need to restart the kernel to use updated packages.\n",
+      "Note: you may need to restart the kernel to use updated packages.\n"
+     ]
+    }
+   ],
+   "source": [
+    "%pip install --upgrade pip ipywidgets -q\n",
+    "%pip install pybop -q\n",
+    "\n",
+    "import time\n",
+    "\n",
+    "import numpy as np\n",
+    "import pybamm\n",
+    "\n",
+    "import pybop"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's fix the random seed in order to generate consistent output during development, although this does not need to be done in practice."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(8)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "id": "5XU-dMtU04p2"
+   },
+   "source": [
+    "### Setting up the model, and problem\n",
+    "\n",
+    "We start by constructing a pybop model, and a synthetic dataset needed for the pybop problem we will be using for the solver benchmarking "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Model\n",
+    "parameter_set = pybop.ParameterSet.pybamm(\"Chen2020\")\n",
+    "model = pybop.lithium_ion.SPM(parameter_set=parameter_set)\n",
+    "\n",
+    "# Synthetic data\n",
+    "t_eval = np.arange(0, 900, 2)\n",
+    "values = model.predict(t_eval=t_eval)\n",
+    "\n",
+    "# Dataset\n",
+    "dataset = pybop.Dataset(\n",
+    "    {\n",
+    "        \"Time [s]\": t_eval,\n",
+    "        \"Current function [A]\": values[\"Current [A]\"].data,\n",
+    "        \"Voltage [V]\": values[\"Voltage [V]\"].data,\n",
+    "    }\n",
+    ")\n",
+    "\n",
+    "# Parameters\n",
+    "parameters = pybop.Parameters(\n",
+    "    pybop.Parameter(\n",
+    "        \"Negative electrode active material volume fraction\",\n",
+    "        prior=pybop.Gaussian(0.6, 0.02),\n",
+    "        bounds=[0.5, 0.8],\n",
+    "    ),\n",
+    "    pybop.Parameter(\n",
+    "        \"Positive electrode active material volume fraction\",\n",
+    "        prior=pybop.Gaussian(0.48, 0.02),\n",
+    "        bounds=[0.4, 0.7],\n",
+    "    ),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "id": "n4OHa-aF04qA"
+   },
+   "source": [
+    "### Defining the solvers for benchmarking\n",
+    "\n",
+    "Now that we have set up the majority of the pybop objects, we construct the solvers we want to benchmark on the given model, and applied current."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "solvers = [\n",
+    "    pybamm.IDAKLUSolver(atol=1e-6, rtol=1e-6),\n",
+    "    pybamm.CasadiSolver(atol=1e-6, rtol=1e-6, mode=\"safe\"),\n",
+    "    pybamm.CasadiSolver(atol=1e-6, rtol=1e-6, mode=\"fast\"),\n",
+    "    pybamm.CasadiSolver(atol=1e-6, rtol=1e-6, mode=\"fast with events\"),\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we construct a range of inputs for the parameters defined above, and select the number of instances in that range to benchmark on. For more statistically repeatable results, increase the variable `n` below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n = 50  # Number of solves\n",
+    "inputs = list(zip(np.linspace(0.45, 0.6, n), np.linspace(0.45, 0.6, n)))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, let's benchmark the solvers without sensitivities. This provides a reference for the non-gradient based pybop optimisers and samplers. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Time Evaluate IDA KLU solver: 0.477\n",
+      "Time Evaluate CasADi solver with 'safe' mode: 0.464\n",
+      "Time Evaluate CasADi solver with 'fast' mode: 0.446\n",
+      "Time Evaluate CasADi solver with 'fast with events' mode: 0.443\n"
+     ]
+    }
+   ],
+   "source": [
+    "for solver in solvers:\n",
+    "    model.solver = solver\n",
+    "    problem = pybop.FittingProblem(model, parameters, dataset)\n",
+    "\n",
+    "    start_time = time.time()\n",
+    "    for input_values in inputs:\n",
+    "        problem.evaluate(inputs=input_values)\n",
+    "    print(f\"Time Evaluate {solver.name}: {time.time() - start_time:.3f}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Excellent, given the above results, we know which solver we should select for optimisation on your machine, i.e. the one with the smallest time. \n",
+    "\n",
+    "Next, let's repeat the same toy problem, but for the gradient-based cost evaluation,"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Time EvaluateS1 IDA KLU solver: 0.997\n",
+      "Time EvaluateS1 CasADi solver with 'safe' mode: 2.743\n",
+      "Time EvaluateS1 CasADi solver with 'fast' mode: 2.689\n",
+      "Time EvaluateS1 CasADi solver with 'fast with events' mode: 2.688\n"
+     ]
+    }
+   ],
+   "source": [
+    "for solver in solvers:\n",
+    "    model.solver = solver\n",
+    "    problem = pybop.FittingProblem(model, parameters, dataset)\n",
+    "\n",
+    "    start_time = time.time()\n",
+    "    for input_values in inputs:\n",
+    "        problem.evaluateS1(inputs=input_values)\n",
+    "    print(f\"Time EvaluateS1 {solver.name}: {time.time() - start_time:.3f}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we have the relevant information for the gradient-based optimisers. Likewise to the above results, we should select the solver with the smallest time."
+   ]
+  }
+ ],
+ "metadata": {
+  "colab": {
+   "provenance": []
+  },
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/notebooks/spm_AdamW.ipynb b/examples/notebooks/spm_AdamW.ipynb
index 9d683d38a..13c61a2cc 100644
--- a/examples/notebooks/spm_AdamW.ipynb
+++ b/examples/notebooks/spm_AdamW.ipynb
@@ -16,7 +16,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/notebooks/spm_electrode_design.ipynb b/examples/notebooks/spm_electrode_design.ipynb
index 6f07b487b..e41cd161a 100644
--- a/examples/notebooks/spm_electrode_design.ipynb
+++ b/examples/notebooks/spm_electrode_design.ipynb
@@ -18,7 +18,7 @@
     "\n",
     "### Setting up the Environment\n",
     "\n",
-    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP from its development branch and upgrade some dependencies:"
+    "Before we begin, we need to ensure that we have all the necessary tools. We will install PyBOP and upgrade dependencies:"
    ]
   },
   {
diff --git a/examples/scripts/selecting_a_solver.py b/examples/scripts/selecting_a_solver.py
new file mode 100644
index 000000000..248d943f5
--- /dev/null
+++ b/examples/scripts/selecting_a_solver.py
@@ -0,0 +1,60 @@
+import time
+
+import numpy as np
+import pybamm
+
+import pybop
+
+# Parameter set and model definition
+parameter_set = pybop.ParameterSet.pybamm("Chen2020")
+model = pybop.lithium_ion.SPM(parameter_set=parameter_set)
+
+solvers = [
+    pybamm.IDAKLUSolver(atol=1e-6, rtol=1e-6),
+    pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-6),
+    pybamm.CasadiSolver(mode="fast with events", atol=1e-6, rtol=1e-6),
+]
+
+# Fitting parameters
+parameters = pybop.Parameters(
+    pybop.Parameter(
+        "Negative electrode active material volume fraction", initial_value=0.55
+    ),
+    pybop.Parameter(
+        "Positive electrode active material volume fraction", initial_value=0.55
+    ),
+)
+
+# Define test protocol and generate data
+experiment = pybop.Experiment([("Discharge at 0.5C for 10 minutes (3 second period)")])
+values = model.predict(
+    initial_state={"Initial open-circuit voltage [V]": 4.2}, experiment=experiment
+)
+
+# Form dataset
+dataset = pybop.Dataset(
+    {
+        "Time [s]": values["Time [s]"].data,
+        "Current function [A]": values["Current [A]"].data,
+        "Voltage [V]": values["Voltage [V]"].data,
+    }
+)
+
+# Create the list of input dicts
+n = 150  # Number of solves
+inputs = list(zip(np.linspace(0.45, 0.6, n), np.linspace(0.45, 0.6, n)))
+
+# Iterate over the solvers and print benchmarks
+for solver in solvers:
+    model.solver = solver
+    problem = pybop.FittingProblem(model, parameters, dataset)
+
+    start_time = time.time()
+    for input_values in inputs:
+        problem.evaluate(inputs=input_values)
+    print(f"Time Evaluate {solver.name}: {time.time() - start_time:.3f}")
+
+    start_time = time.time()
+    for input_values in inputs:
+        problem.evaluateS1(inputs=input_values)
+    print(f"Time EvaluateS1 {solver.name}: {time.time() - start_time:.3f}")
diff --git a/pybop/models/base_model.py b/pybop/models/base_model.py
index 932587cf5..47ae1e076 100644
--- a/pybop/models/base_model.py
+++ b/pybop/models/base_model.py
@@ -792,3 +792,7 @@ def spatial_methods(self):
     @property
     def solver(self):
         return self._solver
+
+    @solver.setter
+    def solver(self, solver):
+        self._solver = solver.copy() if solver is not None else None
diff --git a/pybop/problems/base_problem.py b/pybop/problems/base_problem.py
index 72fc4c197..5b9f4ec11 100644
--- a/pybop/problems/base_problem.py
+++ b/pybop/problems/base_problem.py
@@ -1,5 +1,8 @@
 from typing import Optional
 
+import numpy as np
+from pybamm import IDAKLUSolver
+
 from pybop import BaseModel, Dataset, Parameter, Parameters
 from pybop.parameters.parameter import Inputs
 
@@ -68,6 +71,21 @@ def __init__(
         self._time_data = None
         self._target = None
         self.verbose = False
+        self.failure_output = np.asarray([np.inf])
+
+        # If model.solver is IDAKLU, set output vars for improved performance
+        self.output_vars = tuple(self.signal + self.additional_variables)
+        if self._model is not None and isinstance(self._model.solver, IDAKLUSolver):
+            self._solver_copy = self._model.solver.copy()
+            self._model.solver = IDAKLUSolver(
+                atol=self._solver_copy.atol,
+                rtol=self._solver_copy.rtol,
+                root_method=self._solver_copy.root_method,
+                root_tol=self._solver_copy.root_tol,
+                extrap_tol=self._solver_copy.extrap_tol,
+                options=self._solver_copy._options,  # noqa: SLF001
+                output_variables=self.output_vars,
+            )
 
     def set_initial_state(self, initial_state: Optional[dict] = None):
         """
diff --git a/pybop/problems/fitting_problem.py b/pybop/problems/fitting_problem.py
index 50ca2cf61..252a52d43 100644
--- a/pybop/problems/fitting_problem.py
+++ b/pybop/problems/fitting_problem.py
@@ -131,7 +131,7 @@ def evaluate(
         except Exception as e:
             if self.verbose:
                 print(f"Simulation error: {e}")
-            return {signal: [np.inf] for signal in self.signal}
+            return {signal: self.failure_output for signal in self.signal}
 
         return {
             signal: sol[signal].data
@@ -162,7 +162,9 @@ def evaluateS1(self, inputs: Inputs):
             )
         except Exception as e:
             print(f"Error: {e}")
-            return {signal: [np.inf] for signal in self.signal}, [np.inf]
+            return {
+                signal: self.failure_output for signal in self.signal
+            }, self.failure_output
 
         y = {signal: sol[signal].data for signal in self.signal}
 
diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py
index c441c2b25..2323f0ddc 100644
--- a/tests/integration/test_spm_parameterisations.py
+++ b/tests/integration/test_spm_parameterisations.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pybamm
 import pytest
 
 import pybop
@@ -58,14 +59,29 @@ def init_soc(self, request):
             pybop.MAP,
         ]
     )
-    def cost_class(self, request):
+    def cost(self, request):
         return request.param
 
     def noise(self, sigma, values):
         return np.random.normal(0, sigma, values)
 
+    @pytest.fixture(
+        params=[
+            pybop.SciPyDifferentialEvolution,
+            pybop.AdamW,
+            pybop.CMAES,
+            pybop.CuckooSearch,
+            pybop.IRPropMin,
+            pybop.NelderMead,
+            pybop.SNES,
+            pybop.XNES,
+        ]
+    )
+    def optimiser(self, request):
+        return request.param
+
     @pytest.fixture
-    def spm_costs(self, model, parameters, cost_class, init_soc):
+    def optim(self, optimiser, model, parameters, cost, init_soc):
         # Form dataset
         solution = self.get_data(model, init_soc)
         dataset = pybop.Dataset(
@@ -77,65 +93,64 @@ def spm_costs(self, model, parameters, cost_class, init_soc):
             }
         )
 
-        # Define the cost to optimise
+        # IDAKLU Solver for Gradient-based optimisers
+        if optimiser in [pybop.AdamW, pybop.IRPropMin]:
+            model.solver = pybamm.IDAKLUSolver()
+
+        # Define the problem
         problem = pybop.FittingProblem(model, parameters, dataset)
-        if cost_class in [pybop.GaussianLogLikelihoodKnownSigma]:
-            return cost_class(problem, sigma0=self.sigma0)
-        elif cost_class in [pybop.GaussianLogLikelihood]:
-            return cost_class(problem, sigma0=self.sigma0 * 4)  # Initial sigma0 guess
-        elif cost_class in [pybop.MAP]:
-            return cost_class(
+
+        # Construct the cost
+        if cost in [pybop.GaussianLogLikelihoodKnownSigma]:
+            cost = cost(problem, sigma0=self.sigma0)
+        elif cost in [pybop.GaussianLogLikelihood]:
+            cost = cost(problem, sigma0=self.sigma0 * 4)  # Initial sigma0 guess
+        elif cost in [pybop.MAP]:
+            cost = cost(
                 problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0
             )
-        elif cost_class in [pybop.SumofPower, pybop.Minkowski]:
-            return cost_class(problem, p=2)
+        elif cost in [pybop.SumofPower, pybop.Minkowski]:
+            cost = cost(problem, p=2)
         else:
-            return cost_class(problem)
+            cost = cost(problem)
 
-    @pytest.mark.parametrize(
-        "optimiser",
-        [
-            pybop.SciPyDifferentialEvolution,
-            pybop.AdamW,
-            pybop.CMAES,
-            pybop.CuckooSearch,
-            pybop.IRPropMin,
-            pybop.NelderMead,
-            pybop.SNES,
-            pybop.XNES,
-        ],
-    )
-    @pytest.mark.integration
-    def test_spm_optimisers(self, optimiser, spm_costs):
-        x0 = spm_costs.parameters.initial_value()
+        # Construct optimisation object
         common_args = {
-            "cost": spm_costs,
+            "cost": cost,
             "max_iterations": 250,
             "absolute_tolerance": 1e-6,
             "max_unchanged_iterations": 55,
         }
 
-        # Add sigma0 to ground truth for GaussianLogLikelihood
-        if isinstance(spm_costs, pybop.GaussianLogLikelihood):
-            self.ground_truth = np.concatenate(
-                (self.ground_truth, np.asarray([self.sigma0]))
-            )
-        if isinstance(spm_costs, pybop.MAP):
-            for i in spm_costs.parameters.keys():
-                spm_costs.parameters[i].prior = pybop.Uniform(
+        if isinstance(cost, pybop.MAP):
+            for i in cost.parameters.keys():
+                cost.parameters[i].prior = pybop.Uniform(
                     0.2, 2.0
                 )  # Increase range to avoid prior == np.inf
+
         # Set sigma0 and create optimiser
-        sigma0 = 0.05 if isinstance(spm_costs, pybop.MAP) else None
+        sigma0 = 0.05 if isinstance(cost, pybop.MAP) else None
         optim = optimiser(sigma0=sigma0, **common_args)
 
         # AdamW will use lowest sigma0 for learning rate, so allow more iterations
         if issubclass(optimiser, (pybop.AdamW, pybop.IRPropMin)) and isinstance(
-            spm_costs, pybop.GaussianLogLikelihood
+            cost, pybop.GaussianLogLikelihood
         ):
             common_args["max_unchanged_iterations"] = 75
             optim = optimiser(**common_args)
 
+        return optim
+
+    @pytest.mark.integration
+    def test_spm_optimisers(self, optim):
+        x0 = optim.parameters.initial_value()
+
+        # Add sigma0 to ground truth for GaussianLogLikelihood
+        if isinstance(optim.cost, pybop.GaussianLogLikelihood):
+            self.ground_truth = np.concatenate(
+                (self.ground_truth, np.asarray([self.sigma0]))
+            )
+
         initial_cost = optim.cost(x0)
         x, final_cost = optim.run()
 
@@ -143,7 +158,7 @@ def test_spm_optimisers(self, optimiser, spm_costs):
         if np.allclose(x0, self.ground_truth, atol=1e-5):
             raise AssertionError("Initial guess is too close to ground truth")
 
-        if isinstance(spm_costs, pybop.GaussianLogLikelihood):
+        if isinstance(optim.cost, pybop.GaussianLogLikelihood):
             np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2)
             np.testing.assert_allclose(x[-1], self.sigma0, atol=5e-4)
         else:
@@ -155,7 +170,7 @@ def test_spm_optimisers(self, optimiser, spm_costs):
             np.testing.assert_allclose(x, self.ground_truth, atol=1.5e-2)
 
     @pytest.fixture
-    def spm_two_signal_cost(self, parameters, model, cost_class):
+    def spm_two_signal_cost(self, parameters, model, cost):
         # Form dataset
         solution = self.get_data(model, 0.5)
         dataset = pybop.Dataset(
@@ -175,14 +190,14 @@ def spm_two_signal_cost(self, parameters, model, cost_class):
         signal = ["Voltage [V]", "Bulk open-circuit voltage [V]"]
         problem = pybop.FittingProblem(model, parameters, dataset, signal=signal)
 
-        if cost_class in [pybop.GaussianLogLikelihoodKnownSigma]:
-            return cost_class(problem, sigma0=self.sigma0)
-        elif cost_class in [pybop.MAP]:
-            return cost_class(
+        if cost in [pybop.GaussianLogLikelihoodKnownSigma]:
+            return cost(problem, sigma0=self.sigma0)
+        elif cost in [pybop.MAP]:
+            return cost(
                 problem, pybop.GaussianLogLikelihoodKnownSigma, sigma0=self.sigma0
             )
         else:
-            return cost_class(problem)
+            return cost(problem)
 
     @pytest.mark.parametrize(
         "multi_optimiser",
diff --git a/tests/unit/test_solvers.py b/tests/unit/test_solvers.py
new file mode 100644
index 000000000..3fe318378
--- /dev/null
+++ b/tests/unit/test_solvers.py
@@ -0,0 +1,75 @@
+import numpy as np
+import pybamm
+import pytest
+
+import pybop
+
+
+class TestSolvers:
+    """
+    A class to test the forward model solver interface
+    """
+
+    @pytest.fixture(
+        params=[
+            pybamm.IDAKLUSolver(atol=1e-4, rtol=1e-4),
+            pybamm.CasadiSolver(atol=1e-4, rtol=1e-4, mode="safe"),
+            pybamm.CasadiSolver(atol=1e-4, rtol=1e-4, mode="fast with events"),
+        ]
+    )
+    def solver(self, request):
+        solver = request.param
+        return solver.copy()
+
+    @pytest.fixture
+    def model(self, solver):
+        parameter_set = pybop.ParameterSet.pybamm("Marquis2019")
+        model = pybop.lithium_ion.SPM(parameter_set=parameter_set, solver=solver)
+        return model
+
+    @pytest.mark.unit
+    def test_solvers_with_model_predict(self, model, solver):
+        assert model.solver == solver
+        assert model.solver.atol == 1e-4
+        assert model.solver.rtol == 1e-4
+
+        # Ensure solver is functional
+        sol = model.predict(t_eval=np.linspace(0, 1, 100))
+        assert np.isfinite(sol["Voltage [V]"].data).all()
+
+        signals = ["Voltage [V]", "Bulk open-circuit voltage [V]"]
+        additional_vars = [
+            "Maximum negative particle concentration",
+            "Positive electrode volume-averaged concentration [mol.m-3]",
+        ]
+
+        parameters = pybop.Parameters(
+            pybop.Parameter(
+                "Negative electrode conductivity [S.m-1]", prior=pybop.Uniform(0.1, 100)
+            )
+        )
+        dataset = pybop.Dataset(
+            {
+                "Time [s]": sol["Time [s]"].data,
+                "Current function [A]": sol["Current [A]"].data,
+                "Voltage [V]": sol["Voltage [V]"].data,
+                "Bulk open-circuit voltage [V]": sol[
+                    "Bulk open-circuit voltage [V]"
+                ].data,
+            }
+        )
+        problem = pybop.FittingProblem(
+            model,
+            parameters=parameters,
+            dataset=dataset,
+            signal=signals,
+            additional_variables=additional_vars,
+        )
+
+        y = problem.evaluate(inputs={"Negative electrode conductivity [S.m-1]": 10})
+
+        for signal in signals:
+            assert np.isfinite(y[signal].data).all()
+
+        if isinstance(model.solver, pybamm.IDAKLUSolver):
+            assert model.solver.output_variables is not None