diff --git a/CHANGELOG.md b/CHANGELOG.md
index d18b80e0a8..fc0c452b76 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -294,7 +294,7 @@ example notebook ([#1602](https://github.com/pybamm-team/PyBaMM/pull/1602))
 ## Breaking changes
 
 -   Refactored the `particle` submodel module, with the models having no size distribution now found in `particle.no_distribution`, and those with a size distribution in `particle.size_distribution`. Renamed submodels to indicate the transport model (Fickian diffusion, polynomial profile) and if they are "x-averaged". E.g., `FickianManyParticles` and `FickianSingleParticle` are now `no_distribution.FickianDiffusion` and `no_distribution.XAveragedFickianDiffusion` ([#1602](https://github.com/pybamm-team/PyBaMM/pull/1602))
--   Changed sensitivity API. Removed `ProcessedSymbolicVariable`, all sensitivity now handled within the solvers and `ProcessedVariable` ()  ([#1552](https://github.com/pybamm-team/PyBaMM/pull/1552))
+-   Changed sensitivity API. Removed `ProcessedSymbolicVariable`, all sensitivity now handled within the solvers and `ProcessedVariable` ([#1552](https://github.com/pybamm-team/PyBaMM/pull/1552),[#2276](https://github.com/pybamm-team/PyBaMM/pull/2276))
 -   The `Yang2017` parameter set has been removed as the complete parameter set is not publicly available in the literature ([#1577](https://github.com/pybamm-team/PyBaMM/pull/1577))
 -   Changed how options are specified for the "loss of active material" and "particle cracking" submodels. "loss of active material" can now be one of "none", "stress-driven", or "reaction-driven", or a 2-tuple for different options in negative and positive electrode. Similarly "particle cracking" (now called "particle mechanics") can now be "none", "swelling only", "swelling and cracking", or a 2-tuple ([#1490](https://github.com/pybamm-team/PyBaMM/pull/1490))
 -   Changed the variable in the full diffusion model from "Electrolyte concentration" to "Porosity times concentration" ([#1476](https://github.com/pybamm-team/PyBaMM/pull/1476))
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index c9d05cc53f..1d60418db4 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -200,7 +200,6 @@
 #
 from .solvers.solution import Solution, make_cycle_solution
 from .solvers.processed_variable import ProcessedVariable
-from .solvers.processed_symbolic_variable import ProcessedSymbolicVariable
 from .solvers.base_solver import BaseSolver
 from .solvers.dummy_solver import DummySolver
 from .solvers.algebraic_solver import AlgebraicSolver
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 15374a5ad5..981e22666c 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -1481,20 +1481,11 @@ def _set_up_ext_and_inputs(
         inputs = inputs or {}
 
         # Go through all input parameters that can be found in the model
-        # If any of them are *not* provided by "inputs", a symbolic input parameter is
-        # created, with appropriate size
+        # If any of them are *not* provided by "inputs", raise an error
         for input_param in model.input_parameters:
             name = input_param.name
             if name not in inputs:
-                # Only allow symbolic inputs for CasadiSolver and CasadiAlgebraicSolver
-                if not isinstance(
-                    self, (pybamm.CasadiSolver, pybamm.CasadiAlgebraicSolver)
-                ):
-                    raise pybamm.SolverError(
-                        "Only CasadiSolver and CasadiAlgebraicSolver "
-                        "can have symbolic inputs"
-                    )
-                inputs[name] = casadi.MX.sym(name, input_param._expected_size)
+                raise pybamm.SolverError(f"No value provided for input '{name}'")
 
         external_variables = external_variables or {}
 
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
index 144d369d21..b13c08c27f 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -51,19 +51,10 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         t_eval : :class:`numpy.array`, size (k,)
             The times at which to compute the solution
         inputs_dict : dict, optional
-            Any input parameters to pass to the model when solving. If any input
-            parameters that are present in the model are missing from "inputs", then
-            the solution will consist of `ProcessedSymbolicVariable` objects, which must
-            be provided with inputs to obtain their value.
+            Any input parameters to pass to the model when solving.
         """
         # Record whether there are any symbolic inputs
         inputs_dict = inputs_dict or {}
-        has_symbolic_inputs = any(
-            isinstance(v, casadi.MX) for v in inputs_dict.values()
-        )
-        symbolic_inputs = casadi.vertcat(
-            *[v for v in inputs_dict.values() if isinstance(v, casadi.MX)]
-        )
 
         # Create casadi objects for the root-finder
         inputs = casadi.vertcat(*[v for v in inputs_dict.values()])
@@ -94,7 +85,6 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         y_alg_sym = casadi.MX.sym("y_alg", y0_alg.shape[0])
         y_sym = casadi.vertcat(y0_diff, y_alg_sym)
 
-        t_and_inputs_sym = casadi.vertcat(t_sym, symbolic_inputs)
         alg = model.casadi_algebraic(t_sym, y_sym, inputs)
 
         # Check interpolant extrapolation
@@ -139,7 +129,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         roots = casadi.rootfinder(
             "roots",
             "newton",
-            dict(x=y_alg_sym, p=t_and_inputs_sym, g=alg),
+            dict(x=y_alg_sym, p=t_sym, g=alg),
             {
                 **self.extra_options,
                 "abstol": self.tol,
@@ -150,11 +140,10 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         timer = pybamm.Timer()
         integration_time = 0
         for idx, t in enumerate(t_eval):
-            t_eval_inputs_sym = casadi.vertcat(t, symbolic_inputs)
             # Solve
             try:
                 timer.reset()
-                y_alg_sol = roots(y0_alg, t_eval_inputs_sym)
+                y_alg_sol = roots(y0_alg, t)
                 integration_time += timer.time()
                 success = True
                 message = None
@@ -169,8 +158,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
             # If there are no symbolic inputs, check the function is below the tol
             # Skip this check if there are symbolic inputs
             if success and (
-                has_symbolic_inputs is True
-                or (not any(np.isnan(fun)) and np.all(casadi.fabs(fun) < self.tol))
+                (not any(np.isnan(fun)) and np.all(casadi.fabs(fun) < self.tol))
             ):
                 # update initial guess for the next iteration
                 y0_alg = y_alg_sol
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 7baaa8691c..099af6630c 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -125,19 +125,12 @@ def _integrate(self, model, t_eval, inputs_dict=None):
 
         # Record whether there are any symbolic inputs
         inputs_dict = inputs_dict or {}
-        has_symbolic_inputs = any(
-            isinstance(v, casadi.MX) for v in inputs_dict.values()
-        )
 
         # convert inputs to casadi format
         inputs = casadi.vertcat(*[x for x in inputs_dict.values()])
 
         # Calculate initial event signs needed for some of the modes
-        if (
-            has_symbolic_inputs is False
-            and self.mode != "fast"
-            and model.terminate_events_eval
-        ):
+        if self.mode != "fast" and model.terminate_events_eval:
             init_event_signs = np.sign(
                 np.concatenate(
                     [
@@ -149,18 +142,6 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         else:
             init_event_signs = np.sign([])
 
-        if has_symbolic_inputs:
-            # Create integrator without grid to avoid having to create several times
-            self.create_integrator(model, inputs)
-            solution = self._run_integrator(
-                model,
-                model.y0,
-                inputs_dict,
-                inputs,
-                t_eval,
-                use_grid=False,
-            )
-
         if self.mode in ["fast", "fast with events"] or not model.events:
             if not model.events:
                 pybamm.logger.info("No events found, running fast mode")
diff --git a/pybamm/solvers/processed_symbolic_variable.py b/pybamm/solvers/processed_symbolic_variable.py
deleted file mode 100644
index 76811157c6..0000000000
--- a/pybamm/solvers/processed_symbolic_variable.py
+++ /dev/null
@@ -1,232 +0,0 @@
-#
-# Processed Variable class
-#
-import casadi
-import numbers
-import numpy as np
-
-
-class ProcessedSymbolicVariable(object):
-    """
-    An object that can be evaluated at arbitrary (scalars or vectors) t and x, and
-    returns the (interpolated) value of the base variable at that t and x.
-
-    Parameters
-    ----------
-    base_variable : :class:`pybamm.Symbol`
-        A base variable with a method `evaluate(t,y)` that returns the value of that
-        variable. Note that this can be any kind of node in the expression tree, not
-        just a :class:`pybamm.Variable`.
-        When evaluated, returns an array of size (m,n)
-    solution : :class:`pybamm.Solution`
-        The solution object to be used to create the processed variables
-    """
-
-    def __init__(self, base_variable, solution):
-        # Convert variable to casadi
-        t_MX = casadi.MX.sym("t")
-        y_MX = casadi.MX.sym("y", solution.all_ys[0].shape[0])
-        # Make all inputs symbolic first for converting to casadi
-        all_inputs_as_MX_dict = {}
-        symbolic_inputs_dict = {}
-        for key, value in solution.all_inputs[0].items():
-            if not isinstance(value, casadi.MX):
-                all_inputs_as_MX_dict[key] = casadi.MX.sym("input")
-            else:
-                all_inputs_as_MX_dict[key] = value
-                # Only add symbolic inputs to the "symbolic_inputs" dict
-                symbolic_inputs_dict[key] = value
-
-        all_inputs_as_MX = casadi.vertcat(*[p for p in all_inputs_as_MX_dict.values()])
-        # The symbolic_inputs dictionary will be used for sensitivity
-        symbolic_inputs = casadi.vertcat(*[p for p in symbolic_inputs_dict.values()])
-        var = base_variable.to_casadi(t_MX, y_MX, inputs=all_inputs_as_MX_dict)
-
-        self.base_variable = casadi.Function(
-            "variable", [t_MX, y_MX, all_inputs_as_MX], [var]
-        )
-        # Store some attributes
-        self.t_pts = solution.t
-        self.all_ts = solution.all_ts
-        self.all_ys = solution.all_ys
-        self.mesh = base_variable.mesh
-        self.all_inputs_casadi = solution.all_inputs_casadi
-
-        self.symbolic_inputs_dict = symbolic_inputs_dict
-        self.symbolic_inputs_total_shape = symbolic_inputs.shape[0]
-        self.domain = base_variable.domain
-
-        self.base_eval = self.base_variable(
-            self.all_ts[0][0], self.all_ys[0][:, 0], self.all_inputs_casadi[0]
-        )
-
-        if (
-            isinstance(self.base_eval, numbers.Number)
-            or len(self.base_eval.shape) == 0
-            or self.base_eval.shape[0] == 1
-        ):
-            self.initialise_0D()
-        else:
-            n = self.mesh.npts
-            base_shape = self.base_eval.shape[0]
-            # Try shape that could make the variable a 1D variable
-            if base_shape == n:
-                self.initialise_1D()
-            else:
-                # Raise error for 2D variable
-                raise NotImplementedError(
-                    "Shape not recognized for {} ".format(base_variable)
-                    + "(note processing of 2D and 3D variables is not yet "
-                    + "implemented)"
-                )
-
-        # Make entries a function and compute jacobian
-        entries_MX = self.entries
-        self.casadi_entries_fn = casadi.Function(
-            "variable", [symbolic_inputs], [entries_MX]
-        )
-
-        # Don't compute jacobian if the entries are a DM (not symbolic)
-        if isinstance(entries_MX, casadi.DM):
-            self.casadi_sens_fn = None
-        # Do compute jacobian if the entries are symbolic (functions of input)
-        else:
-            sens_MX = casadi.jacobian(entries_MX, symbolic_inputs)
-            self.casadi_sens_fn = casadi.Function(
-                "variable", [symbolic_inputs], [sens_MX]
-            )
-
-    def initialise_0D(self):
-        """Create a 0D variable"""
-        # Evaluate the base_variable index-by-index
-        idx = 0
-        for ts, ys, inputs in zip(self.all_ts, self.all_ys, self.all_inputs_casadi):
-            for inner_idx, t in enumerate(ts):
-                t = ts[inner_idx]
-                y = ys[:, inner_idx]
-                next_entries = self.base_variable(t, y, inputs)
-                if idx == 0:
-                    entries = next_entries
-                else:
-                    entries = casadi.horzcat(entries, next_entries)
-                idx += 1
-
-        self.entries = entries
-        self.dimensions = 0
-
-    def initialise_1D(self):
-        """Create a 1D variable"""
-        len_space = self.base_eval.shape[0]
-        entries = np.empty((len_space, len(self.t_pts)))
-
-        # Evaluate the base_variable index-by-index
-        idx = 0
-        for ts, ys, inputs in zip(self.all_ts, self.all_ys, self.all_inputs_casadi):
-            for inner_idx, t in enumerate(ts):
-                t = ts[inner_idx]
-                y = ys[:, inner_idx]
-                next_entries = self.base_variable(t, y, inputs)
-                if idx == 0:
-                    entries = next_entries
-                else:
-                    entries = casadi.vertcat(entries, next_entries)
-                idx += 1
-
-        self.entries = entries
-
-    def value(self, inputs=None, check_inputs=True):
-        """
-        Returns the value of the variable at the specified input values
-
-        Parameters
-        ----------
-        inputs : dict
-            The inputs at which to evaluate the variable.
-
-        Returns
-        -------
-        casadi.DM
-            A casadi matrix of size (n_x * n_t, 1), where n_x is the number of spatial
-            discretisation points for the variable, and n_t is the length of the time
-            vector
-        """
-        if inputs is None:
-            return self.casadi_entries_fn(casadi.DM())
-        else:
-            if check_inputs:
-                inputs = self._check_and_transform(inputs)
-            return self.casadi_entries_fn(inputs)
-
-    def sensitivity(self, inputs=None, check_inputs=True):
-        """
-        Returns the sensitivity of the variable to the symbolic inputs at the specified
-        input values
-
-        Parameters
-        ----------
-        inputs : dict
-            The inputs at which to evaluate the variable.
-
-        Returns
-        -------
-        casadi.DM
-            A casadi matrix of size (n_x * n_t, n_p), where n_x is the number of spatial
-            discretisation points for the variable, n_t is the length of the time
-            vector, and n_p is the number of input parameters
-        """
-        if self.casadi_sens_fn is None:
-            raise ValueError(
-                "Variable is not symbolic, so sensitivities are not defined"
-            )
-        if check_inputs:
-            inputs = self._check_and_transform(inputs)
-        return self.casadi_sens_fn(inputs)
-
-    def value_and_sensitivity(self, inputs=None):
-        """
-        Returns the value of the variable and its sensitivity to the symbolic inputs at
-        the specified input values
-
-        Parameters
-        ----------
-        inputs : dict
-            The inputs at which to evaluate the variable.
-        """
-        inputs = self._check_and_transform(inputs)
-        # Pass check_inputs=False to avoid re-checking inputs
-        return (
-            self.value(inputs, check_inputs=False),
-            self.sensitivity(inputs, check_inputs=False),
-        )
-
-    def _check_and_transform(self, inputs_dict):
-        """Check dictionary has the right inputs, and convert to a vector"""
-        # Convert dict to casadi vector
-        if not isinstance(inputs_dict, dict):
-            raise TypeError("inputs should be 'dict' but are {}".format(inputs_dict))
-        # Sort input dictionary keys according to the symbolic inputs dictionary
-        # For practical number of input parameters this should be extremely fast and
-        # so is ok to do at each step
-        try:
-            inputs_dict_sorted = {
-                k: inputs_dict[k] for k in self.symbolic_inputs_dict.keys()
-            }
-        except KeyError as e:
-            raise KeyError("Inconsistent input keys. '{}' not found".format(e.args[0]))
-        inputs = casadi.vertcat(*[p for p in inputs_dict_sorted.values()])
-        if inputs.shape[0] != self.symbolic_inputs_total_shape:
-            # Find the variable which caused the error, for a clearer error message
-            for key, inp in inputs_dict_sorted.items():
-                if inp.shape[0] != self.symbolic_inputs_dict[key].shape[0]:
-                    raise ValueError(
-                        "Wrong shape for input '{}': expected {}, actual {}".format(
-                            key, self.symbolic_inputs_dict[key].shape[0], inp.shape[0]
-                        )
-                    )
-
-        return inputs
-
-    @property
-    def data(self):
-        """Same as entries, but different name"""
-        return self.entries
diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py
index e03754cb5f..968e95e1a6 100644
--- a/pybamm/solvers/processed_variable.py
+++ b/pybamm/solvers/processed_variable.py
@@ -47,8 +47,6 @@ def __init__(self, base_variables, base_variables_casadi, solution, warn=True):
         self.domains = base_variables[0].domains
         self.warn = warn
 
-        self.symbolic_inputs = solution.has_symbolic_inputs
-
         # Sensitivity starts off uninitialized, only set when called
         self._sensitivities = None
         self.solution_sensitivities = solution.sensitivities
@@ -307,15 +305,11 @@ def initialise_2D(self):
             self.second_dimension = "x"
             self.r_sol = first_dim_pts
             self.x_sol = second_dim_pts
-        elif (
-            self.domain[0]
-            in [
-                "negative electrode",
-                "separator",
-                "positive electrode",
-            ]
-            and self.domains["secondary"] == ["current collector"]
-        ):
+        elif self.domain[0] in [
+            "negative electrode",
+            "separator",
+            "positive electrode",
+        ] and self.domains["secondary"] == ["current collector"]:
             self.first_dimension = "x"
             self.second_dimension = "z"
             self.x_sol = first_dim_pts
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 992259123e..cbcf74ff06 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -102,12 +102,8 @@ def __init__(
         self._y_event = y_event
         self._termination = termination
 
-        self.has_symbolic_inputs = any(
-            isinstance(v, casadi.MX) for v in self.all_inputs[0].values()
-        )
-
         # Check no ys are too large
-        if check_solution and not self.has_symbolic_inputs:
+        if check_solution:
             self.check_ys_are_not_too_large()
 
         # Copy the timescale_eval and lengthscale_evals if they exist
@@ -476,51 +472,36 @@ def update(self, variables):
         # Process
         for key in variables:
             pybamm.logger.debug("Post-processing {}".format(key))
-            # If there are symbolic inputs then we need to make a
-            # ProcessedSymbolicVariable
-            if self.has_symbolic_inputs is True:
-                var = pybamm.ProcessedSymbolicVariable(
-                    self.all_models[0].variables_and_events[key], self
-                )
-
-            # Otherwise a standard ProcessedVariable is ok
-            else:
-                vars_pybamm = [
-                    model.variables_and_events[key] for model in self.all_models
-                ]
-
-                # Iterate through all models, some may be in the list several times and
-                # therefore only get set up once
-                vars_casadi = []
-                for model, ys, inputs, var_pybamm in zip(
-                    self.all_models, self.all_ys, self.all_inputs, vars_pybamm
-                ):
-                    if key in model._variables_casadi:
-                        var_casadi = model._variables_casadi[key]
-                    else:
-                        t_MX = casadi.MX.sym("t")
-                        y_MX = casadi.MX.sym("y", ys.shape[0])
-                        symbolic_inputs_dict = {
-                            key: casadi.MX.sym("input", value.shape[0])
-                            for key, value in inputs.items()
-                        }
-                        symbolic_inputs = casadi.vertcat(
-                            *[p for p in symbolic_inputs_dict.values()]
-                        )
-
-                        # Convert variable to casadi
-                        # Make all inputs symbolic first for converting to casadi
-                        var_sym = var_pybamm.to_casadi(
-                            t_MX, y_MX, inputs=symbolic_inputs_dict
-                        )
-
-                        var_casadi = casadi.Function(
-                            "variable", [t_MX, y_MX, symbolic_inputs], [var_sym]
-                        )
-                        model._variables_casadi[key] = var_casadi
-                    vars_casadi.append(var_casadi)
+            vars_pybamm = [model.variables_and_events[key] for model in self.all_models]
+
+            # Iterate through all models, some may be in the list several times and
+            # therefore only get set up once
+            vars_casadi = []
+            for model, ys, inputs, var_pybamm in zip(
+                self.all_models, self.all_ys, self.all_inputs, vars_pybamm
+            ):
+                if key in model._variables_casadi:
+                    var_casadi = model._variables_casadi[key]
+                else:
+                    t_MX = casadi.MX.sym("t")
+                    y_MX = casadi.MX.sym("y", ys.shape[0])
+                    inputs_MX_dict = {
+                        key: casadi.MX.sym("input", value.shape[0])
+                        for key, value in inputs.items()
+                    }
+                    inputs_MX = casadi.vertcat(*[p for p in inputs_MX_dict.values()])
+
+                    # Convert variable to casadi
+                    # Make all inputs symbolic first for converting to casadi
+                    var_sym = var_pybamm.to_casadi(t_MX, y_MX, inputs=inputs_MX_dict)
+
+                    var_casadi = casadi.Function(
+                        "variable", [t_MX, y_MX, inputs_MX], [var_sym]
+                    )
+                    model._variables_casadi[key] = var_casadi
+                vars_casadi.append(var_casadi)
 
-                var = pybamm.ProcessedVariable(vars_pybamm, vars_casadi, self)
+            var = pybamm.ProcessedVariable(vars_pybamm, vars_casadi, self)
 
             # Save variable and data
             self._variables[key] = var
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index d05ebd1698..3f184ef114 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -91,7 +91,7 @@ def test_block_symbolic_inputs(self):
         model.rhs = {a: a * p}
         with self.assertRaisesRegex(
             pybamm.SolverError,
-            "Only CasadiSolver and CasadiAlgebraicSolver can have symbolic inputs",
+            "No value provided for input 'p'"
         ):
             solver.solve(model, np.array([1, 2, 3]))
 
@@ -217,7 +217,7 @@ def rhs_eval(self, t, y, inputs):
 
             def algebraic_eval(self, t, y, inputs):
                 # algebraic equation has no root
-                return y ** 2 + 1
+                return y**2 + 1
 
         solver = pybamm.BaseSolver(root_method="hybr")
 
@@ -349,7 +349,7 @@ def exact_diff_b(y, a, b):
             u = pybamm.Variable("u")
             a = pybamm.InputParameter("a")
             b = pybamm.InputParameter("b")
-            model.rhs = {v: a * v ** 2 + b * v + a ** 2}
+            model.rhs = {v: a * v**2 + b * v + a**2}
             model.algebraic = {u: a * v - u}
             model.initial_conditions = {v: 1, u: a * 1}
             model.convert_to_format = convert_to_format
diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index f7af342562..6d432fed9d 100644
--- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -58,13 +58,13 @@ class Model:
             p = casadi.MX.sym("p")
             length_scales = {}
             rhs = {}
-            casadi_algebraic = casadi.Function("alg", [t, y, p], [y ** 2 + 1])
+            casadi_algebraic = casadi.Function("alg", [t, y, p], [y**2 + 1])
             bounds = (np.array([-np.inf]), np.array([np.inf]))
             interpolant_extrapolation_events_eval = []
 
             def algebraic_eval(self, t, y, inputs):
                 # algebraic equation has no real root
-                return y ** 2 + 1
+                return y**2 + 1
 
         model = Model()
 
@@ -86,13 +86,13 @@ class NaNModel:
             y = casadi.MX.sym("y")
             p = casadi.MX.sym("p")
             rhs = {}
-            casadi_algebraic = casadi.Function("alg", [t, y, p], [y ** 0.5])
+            casadi_algebraic = casadi.Function("alg", [t, y, p], [y**0.5])
             bounds = (np.array([-np.inf]), np.array([np.inf]))
             interpolant_extrapolation_events_eval = []
 
             def algebraic_eval(self, t, y, inputs):
                 # algebraic equation has no real root
-                return y ** 0.5
+                return y**0.5
 
         model = NaNModel()
         with self.assertRaisesRegex(
@@ -189,15 +189,22 @@ def test_solve_with_symbolic_input(self):
 
         # Solve
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
-        np.testing.assert_array_equal(solution["var"].value({"param": 7}), -7)
-        np.testing.assert_array_equal(solution["var"].value({"param": 3}), -3)
-        np.testing.assert_array_equal(solution["var"].sensitivity({"param": 3}), -1)
+        solution = solver.solve(
+            model, [0], inputs={"param": 7}, calculate_sensitivities=True
+        )
+        np.testing.assert_array_equal(solution["var"].data, -7)
+
+        solution = solver.solve(
+            model, [0], inputs={"param": 3}, calculate_sensitivities=True
+        )
+        np.testing.assert_array_equal(solution["var"].data, -3)
+        np.testing.assert_array_equal(solution["var"].sensitivities["param"], -1)
 
     def test_least_squares_fit(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var", domain="negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         p = pybamm.InputParameter("p")
         q = pybamm.InputParameter("q")
         model.algebraic = {var: (var - p)}
@@ -210,18 +217,25 @@ def test_least_squares_fit(self):
 
         # Solve
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
-        sol_var = solution["objective"]
 
         def objective(x):
-            return sol_var.value({"p": x[0], "q": x[1]}).full().flatten()
+            solution = solver.solve(
+                model, [0], inputs={"p": x[0], "q": x[1]}, calculate_sensitivities=True
+            )
+            return solution["objective"].data.flatten()
 
         # without jacobian
         lsq_sol = least_squares(objective, [2, 2], method="lm")
         np.testing.assert_array_almost_equal(lsq_sol.x, [3, 3], decimal=3)
 
         def jac(x):
-            return sol_var.sensitivity({"p": x[0], "q": x[1]})
+            solution = solver.solve(
+                model, [0], inputs={"p": x[0], "q": x[1]}, calculate_sensitivities=True
+            )
+            return np.concatenate(
+                [solution["objective"].sensitivities[name] for name in ["p", "q"]],
+                axis=1,
+            )
 
         # with jacobian
         lsq_sol = least_squares(objective, [2, 2], jac=jac, method="lm")
@@ -230,6 +244,7 @@ def jac(x):
     def test_solve_with_symbolic_input_1D_scalar_input(self):
         var = pybamm.Variable("var", "negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         param = pybamm.InputParameter("param")
         model.algebraic = {var: var + param}
         model.initial_conditions = {var: 2}
@@ -241,14 +256,21 @@ def test_solve_with_symbolic_input_1D_scalar_input(self):
 
         # Solve - scalar input
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
-        np.testing.assert_array_equal(solution["var"].value({"param": 7}), -7)
-        np.testing.assert_array_equal(solution["var"].value({"param": 3}), -3)
-        np.testing.assert_array_equal(solution["var"].sensitivity({"param": 3}), -1)
+        solution = solver.solve(
+            model, [0], inputs={"param": 7}, calculate_sensitivities=True
+        )
+        np.testing.assert_array_equal(solution["var"].data, -7)
+
+        solution = solver.solve(
+            model, [0], inputs={"param": 3}, calculate_sensitivities=True
+        )
+        np.testing.assert_array_equal(solution["var"].data, -3)
+        np.testing.assert_array_equal(solution["var"].sensitivities["param"], -1)
 
     def test_solve_with_symbolic_input_1D_vector_input(self):
         var = pybamm.Variable("var", "negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         param = pybamm.InputParameter("param", "negative electrode")
         model.algebraic = {var: var + param}
         model.initial_conditions = {var: 2}
@@ -260,23 +282,23 @@ def test_solve_with_symbolic_input_1D_vector_input(self):
 
         # Solve - scalar input
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
         n = disc.mesh["negative electrode"].npts
 
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
-        p = np.linspace(0, 1, n)[:, np.newaxis]
-        np.testing.assert_array_almost_equal(
-            solution["var"].value({"param": 3 * np.ones(n)}), -3
+        solution = solver.solve(
+            model, [0], inputs={"param": 3 * np.ones(n)}, calculate_sensitivities=True
         )
         np.testing.assert_array_almost_equal(
-            solution["var"].value({"param": 2 * p}), -2 * p
+            solution["var"].sensitivities["param"], -np.eye(40)
         )
-        np.testing.assert_array_almost_equal(
-            solution["var"].sensitivity({"param": 3 * np.ones(n)}), -np.eye(40)
+        np.testing.assert_array_almost_equal(solution["var"].data, -3)
+        p = np.linspace(0, 1, n)[:, np.newaxis]
+        solution = solver.solve(
+            model, [0], inputs={"param": p}, calculate_sensitivities=True
         )
+        np.testing.assert_array_almost_equal(solution["var"].data, -p)
         np.testing.assert_array_almost_equal(
-            solution["var"].sensitivity({"param": p}), -np.eye(40)
+            solution["var"].sensitivities["param"], -np.eye(40)
         )
 
     def test_solve_with_symbolic_input_in_initial_conditions(self):
@@ -293,15 +315,21 @@ def test_solve_with_symbolic_input_in_initial_conditions(self):
 
         # Solve
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
-        np.testing.assert_array_equal(solution["var"].value({"param": 7}), -2)
-        np.testing.assert_array_equal(solution["var"].value({"param": 3}), -2)
-        np.testing.assert_array_equal(solution["var"].sensitivity({"param": 3}), 0)
+        solution = solver.solve(
+            model, [0], inputs={"param": 7}, calculate_sensitivities=True
+        )
+        np.testing.assert_array_equal(solution["var"].data, -2)
+        np.testing.assert_array_equal(solution["var"].sensitivities["param"], 0)
+        solution = solver.solve(
+            model, [0], inputs={"param": 3}, calculate_sensitivities=True
+        )
+        np.testing.assert_array_equal(solution["var"].data, -2)
 
     def test_least_squares_fit_input_in_initial_conditions(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var", domain="negative electrode")
         model = pybamm.BaseModel()
+        model._length_scales = {"negative electrode": pybamm.Scalar(1)}
         p = pybamm.InputParameter("p")
         q = pybamm.InputParameter("q")
         model.algebraic = {var: (var - p)}
@@ -314,11 +342,12 @@ def test_least_squares_fit_input_in_initial_conditions(self):
 
         # Solve
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, [0])
-        sol_var = solution["objective"]
 
         def objective(x):
-            return sol_var.value({"p": x[0], "q": x[1]}).full().flatten()
+            solution = solver.solve(
+                model, [0], inputs={"p": x[0], "q": x[1]}, calculate_sensitivities=True
+            )
+            return solution["objective"].data.flatten()
 
         # without jacobian
         lsq_sol = least_squares(objective, [2, 2], method="lm")
diff --git a/tests/unit/test_solvers/test_processed_symbolic_variable.py b/tests/unit/test_solvers/test_processed_symbolic_variable.py
deleted file mode 100644
index 5f6ee9921f..0000000000
--- a/tests/unit/test_solvers/test_processed_symbolic_variable.py
+++ /dev/null
@@ -1,327 +0,0 @@
-#
-# Tests for the Processed Variable class
-#
-import pybamm
-import casadi
-
-import numpy as np
-import unittest
-import tests
-
-
-class TestProcessedSymbolicVariable(unittest.TestCase):
-    def test_processed_variable_0D(self):
-        # without inputs
-        y = pybamm.StateVector(slice(0, 1))
-        var = 2 * y
-        var.mesh = None
-
-        t_sol = np.linspace(0, 1)
-        y_sol = np.array([np.linspace(0, 5)])
-        solution = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        processed_var = pybamm.ProcessedSymbolicVariable(var, solution)
-        np.testing.assert_array_equal(processed_var.value(), 2 * y_sol)
-
-        # No sensitivity as variable is not symbolic
-        with self.assertRaisesRegex(ValueError, "Variable is not symbolic"):
-            processed_var.sensitivity()
-
-    def test_processed_variable_0D_with_inputs(self):
-        # with symbolic inputs
-        y = pybamm.StateVector(slice(0, 1))
-        p = pybamm.InputParameter("p")
-        q = pybamm.InputParameter("q")
-        var = p * y + q
-        var.mesh = None
-
-        t_sol = np.linspace(0, 1)
-        y_sol = np.array([np.linspace(0, 5)])
-        solution = pybamm.Solution(
-            t_sol,
-            y_sol,
-            pybamm.BaseModel(),
-            {"p": casadi.MX.sym("p"), "q": casadi.MX.sym("q")},
-        )
-        processed_var = pybamm.ProcessedSymbolicVariable(var, solution)
-        np.testing.assert_array_equal(
-            processed_var.value({"p": 3, "q": 4}).full(), 3 * y_sol + 4
-        )
-        np.testing.assert_array_equal(
-            processed_var.sensitivity({"p": 3, "q": 4}).full(),
-            np.c_[y_sol.T, np.ones_like(y_sol).T],
-        )
-
-        # via value_and_sensitivity
-        val, sens = processed_var.value_and_sensitivity({"p": 3, "q": 4})
-        np.testing.assert_array_equal(val.full(), 3 * y_sol + 4)
-        np.testing.assert_array_equal(
-            sens.full(), np.c_[y_sol.T, np.ones_like(y_sol).T]
-        )
-
-        # Test bad inputs
-        with self.assertRaisesRegex(TypeError, "inputs should be 'dict'"):
-            processed_var.value(1)
-        with self.assertRaisesRegex(KeyError, "Inconsistent input keys"):
-            processed_var.value({"not p": 3})
-
-    def test_processed_variable_0D_some_inputs(self):
-        # with some symbolic inputs and some non-symbolic inputs
-        y = pybamm.StateVector(slice(0, 1))
-        p = pybamm.InputParameter("p")
-        q = pybamm.InputParameter("q")
-        var = p * y - q
-        var.mesh = None
-
-        t_sol = np.linspace(0, 1)
-        y_sol = np.array([np.linspace(0, 5)])
-        solution = pybamm.Solution(
-            t_sol, y_sol, pybamm.BaseModel(), {"p": casadi.MX.sym("p"), "q": 2}
-        )
-        processed_var = pybamm.ProcessedSymbolicVariable(var, solution)
-        np.testing.assert_array_equal(
-            processed_var.value({"p": 3}).full(), 3 * y_sol - 2
-        )
-        np.testing.assert_array_equal(
-            processed_var.sensitivity({"p": 3}).full(), y_sol.T
-        )
-
-    def test_processed_variable_1D(self):
-        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
-        x = pybamm.SpatialVariable("x", domain=["negative electrode", "separator"])
-        eqn = var + x
-
-        # On nodes
-        disc = tests.get_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        x_sol = disc.process_symbol(x).entries[:, 0]
-        eqn_sol = disc.process_symbol(eqn)
-
-        # With scalar t_sol
-        t_sol = np.array([0])
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        processed_eqn = pybamm.ProcessedSymbolicVariable(eqn_sol, sol)
-        np.testing.assert_array_equal(
-            processed_eqn.value(), y_sol + x_sol[:, np.newaxis]
-        )
-
-        # With vector t_sol
-        t_sol = np.linspace(0, 1)
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * np.linspace(0, 5)
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        processed_eqn = pybamm.ProcessedSymbolicVariable(eqn_sol, sol)
-        np.testing.assert_array_equal(
-            processed_eqn.value(), (y_sol + x_sol[:, np.newaxis]).T.reshape(-1, 1)
-        )
-
-    def test_processed_variable_1D_with_scalar_inputs(self):
-        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
-        x = pybamm.SpatialVariable("x", domain=["negative electrode", "separator"])
-        p = pybamm.InputParameter("p")
-        q = pybamm.InputParameter("q")
-        eqn = var * p + 2 * q
-
-        # On nodes
-        disc = tests.get_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        x_sol = disc.process_symbol(x).entries[:, 0]
-        eqn_sol = disc.process_symbol(eqn)
-
-        # Scalar t
-        t_sol = np.array([0])
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
-
-        sol = pybamm.Solution(
-            t_sol,
-            y_sol,
-            pybamm.BaseModel(),
-            {"p": casadi.MX.sym("p"), "q": casadi.MX.sym("q")},
-        )
-        processed_eqn = pybamm.ProcessedSymbolicVariable(eqn_sol, sol)
-
-        # Test values
-        np.testing.assert_array_equal(
-            processed_eqn.value({"p": 27, "q": -42}), 27 * y_sol - 84
-        )
-
-        # Test sensitivities
-        np.testing.assert_array_equal(
-            processed_eqn.sensitivity({"p": 27, "q": -84}),
-            np.c_[y_sol, 2 * np.ones_like(y_sol)],
-        )
-
-        ################################################################################
-        # Vector t
-        t_sol = np.linspace(0, 1)
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * np.linspace(0, 5)
-
-        sol = pybamm.Solution(
-            t_sol,
-            y_sol,
-            pybamm.BaseModel(),
-            {"p": casadi.MX.sym("p"), "q": casadi.MX.sym("q")},
-        )
-        processed_eqn = pybamm.ProcessedSymbolicVariable(eqn_sol, sol)
-
-        # Test values
-        np.testing.assert_array_equal(
-            processed_eqn.value({"p": 27, "q": -42}), (27 * y_sol - 84).T.reshape(-1, 1)
-        )
-
-        # Test sensitivities
-        np.testing.assert_array_equal(
-            processed_eqn.sensitivity({"p": 27, "q": -42}),
-            np.c_[y_sol.T.flatten(), 2 * np.ones_like(y_sol.T.flatten())],
-        )
-
-    def test_processed_variable_1D_with_vector_inputs(self):
-        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
-        x = pybamm.SpatialVariable("x", domain=["negative electrode", "separator"])
-        p = pybamm.InputParameter("p", domain=["negative electrode", "separator"])
-        q = pybamm.InputParameter("q")
-        eqn = (var * p) ** 2 + 2 * q
-
-        # On nodes
-        disc = tests.get_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        x_sol = disc.process_symbol(x).entries[:, 0]
-        n = x_sol.size
-        eqn_sol = disc.process_symbol(eqn)
-
-        # Scalar t
-        t_sol = np.array([0])
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(
-            t_sol,
-            y_sol,
-            pybamm.BaseModel(),
-            {"p": casadi.MX.sym("p", n), "q": casadi.MX.sym("q")},
-        )
-        processed_eqn = pybamm.ProcessedSymbolicVariable(eqn_sol, sol)
-
-        # Test values - constant p
-        np.testing.assert_array_equal(
-            processed_eqn.value({"p": 27 * np.ones(n), "q": -42}),
-            (27 * y_sol) ** 2 - 84,
-        )
-        # Test values - varying p
-        p = np.linspace(0, 1, n)
-        np.testing.assert_array_equal(
-            processed_eqn.value({"p": p, "q": 3}), (p[:, np.newaxis] * y_sol) ** 2 + 6
-        )
-
-        # Test sensitivities - constant p
-        np.testing.assert_array_equal(
-            processed_eqn.sensitivity({"p": 2 * np.ones(n), "q": -84}),
-            np.c_[100 * np.eye(y_sol.size), 2 * np.ones(n)],
-        )
-        # Test sensitivities - varying p
-        # d/dy((py)**2) = (2*p*y) * y
-        np.testing.assert_array_equal(
-            processed_eqn.sensitivity({"p": p, "q": -84}),
-            np.c_[
-                np.diag((2 * p[:, np.newaxis] * y_sol ** 2).flatten()), 2 * np.ones(n)
-            ],
-        )
-
-        # Bad shape
-        with self.assertRaisesRegex(
-            ValueError, "Wrong shape for input 'p': expected 65, actual 5"
-        ):
-            processed_eqn.value({"p": casadi.MX.sym("p", 5), "q": 1})
-
-    def test_1D_different_domains(self):
-        # Negative electrode domain
-        var = pybamm.Variable("var", domain=["negative electrode"])
-        x = pybamm.SpatialVariable("x", domain=["negative electrode"])
-
-        disc = tests.get_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        x_sol = disc.process_symbol(x).entries[:, 0]
-        var_sol = disc.process_symbol(var)
-
-        t_sol = np.array([0])
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        pybamm.ProcessedSymbolicVariable(var_sol, sol)
-
-        # Particle domain
-        var = pybamm.Variable("var", domain=["negative particle"])
-        r = pybamm.SpatialVariable("r", domain=["negative particle"])
-
-        disc = tests.get_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        r_sol = disc.process_symbol(r).entries[:, 0]
-        var_sol = disc.process_symbol(var)
-
-        t_sol = np.array([0])
-        y_sol = np.ones_like(r_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        pybamm.ProcessedSymbolicVariable(var_sol, sol)
-
-        # Current collector domain
-        var = pybamm.Variable("var", domain=["current collector"])
-        z = pybamm.SpatialVariable("z", domain=["current collector"])
-
-        disc = tests.get_1p1d_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        z_sol = disc.process_symbol(z).entries[:, 0]
-        var_sol = disc.process_symbol(var)
-
-        t_sol = np.array([0])
-        y_sol = np.ones_like(z_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        pybamm.ProcessedSymbolicVariable(var_sol, sol)
-
-        # Other domain
-        var = pybamm.Variable("var", domain=["line"])
-        x = pybamm.SpatialVariable("x", domain=["line"])
-
-        geometry = pybamm.Geometry(
-            {"line": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
-        )
-        submesh_types = {"line": pybamm.Uniform1DSubMesh}
-        var_pts = {x: 10}
-        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-        disc = pybamm.Discretisation(mesh, {"line": pybamm.FiniteVolume()})
-        disc.set_variable_slices([var])
-        x_sol = disc.process_symbol(x).entries[:, 0]
-        var_sol = disc.process_symbol(var)
-
-        t_sol = np.array([0])
-        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        pybamm.ProcessedSymbolicVariable(var_sol, sol)
-
-        # 2D fails
-        var = pybamm.Variable(
-            "var",
-            domain=["negative particle"],
-            auxiliary_domains={"secondary": "negative electrode"},
-        )
-        r = pybamm.SpatialVariable(
-            "r",
-            domain=["negative particle"],
-            auxiliary_domains={"secondary": "negative electrode"},
-        )
-
-        disc = tests.get_p2d_discretisation_for_testing()
-        disc.set_variable_slices([var])
-        r_sol = disc.process_symbol(r).entries[:, 0]
-        var_sol = disc.process_symbol(var)
-
-        t_sol = np.array([0])
-        y_sol = np.ones_like(r_sol)[:, np.newaxis] * 5
-        sol = pybamm.Solution(t_sol, y_sol, pybamm.BaseModel(), {})
-        with self.assertRaisesRegex(NotImplementedError, "Shape not recognized"):
-            pybamm.ProcessedSymbolicVariable(var_sol, sol)
-
-
-if __name__ == "__main__":
-    print("Add -v for more debug output")
-    import sys
-
-    if "-v" in sys.argv:
-        debug = True
-    pybamm.settings.debug_mode = True
-    unittest.main()
diff --git a/tests/unit/test_solvers/test_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py
index b6814668d6..f0ac25c683 100644
--- a/tests/unit/test_solvers/test_processed_variable.py
+++ b/tests/unit/test_solvers/test_processed_variable.py
@@ -13,16 +13,16 @@ def to_casadi(var_pybamm, y, inputs=None):
     t_MX = casadi.MX.sym("t")
     y_MX = casadi.MX.sym("y", y.shape[0])
 
-    symbolic_inputs_dict = {}
+    inputs_MX_dict = {}
     inputs = inputs or {}
     for key, value in inputs.items():
-        symbolic_inputs_dict[key] = casadi.MX.sym("input", value.shape[0])
+        inputs_MX_dict[key] = casadi.MX.sym("input", value.shape[0])
 
-    symbolic_inputs = casadi.vertcat(*[p for p in symbolic_inputs_dict.values()])
+    inputs_MX = casadi.vertcat(*[p for p in inputs_MX_dict.values()])
 
-    var_sym = var_pybamm.to_casadi(t_MX, y_MX, inputs=symbolic_inputs_dict)
+    var_sym = var_pybamm.to_casadi(t_MX, y_MX, inputs=inputs_MX_dict)
 
-    var_casadi = casadi.Function("variable", [t_MX, y_MX, symbolic_inputs], [var_sym])
+    var_casadi = casadi.Function("variable", [t_MX, y_MX, inputs_MX], [var_sym])
     return var_casadi