diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8638d85ec8..f940e63a82 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,8 @@
 ## Features
 
 -   Added functionality to solver to automatically discretise a 0D model ([#947](https://github.com/pybamm-team/PyBaMM/pull/947))
+-   Added sensitivity to `CasadiAlgebraicSolver` ([#940](https://github.com/pybamm-team/PyBaMM/pull/940))
+-   Added `ProcessedSymbolicVariable` class, which can handle symbolic variables (i.e. variables for which the inputs are symbolic) ([#940](https://github.com/pybamm-team/PyBaMM/pull/940))
 -   Made `QuickPlot` compatible with Google Colab ([#935](https://github.com/pybamm-team/PyBaMM/pull/935))
 -   Added `BasicFull` model for lead-acid ([#932](https://github.com/pybamm-team/PyBaMM/pull/932))
 
diff --git a/docs/index.rst b/docs/index.rst
index e7910fdc97..1ec6349371 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -33,7 +33,6 @@ Contents
     source/experiments/index
     source/simulation
     source/quick_plot
-    source/processed_variable
     source/util
     source/citations
     source/parameters_cli
diff --git a/docs/source/solvers/index.rst b/docs/source/solvers/index.rst
index 7954aea119..1aff323a6a 100644
--- a/docs/source/solvers/index.rst
+++ b/docs/source/solvers/index.rst
@@ -10,3 +10,5 @@ Solvers
   casadi_solver
   algebraic_solvers
   solution
+  processed_variable
+
diff --git a/docs/source/processed_variable.rst b/docs/source/solvers/processed_variable.rst
similarity index 61%
rename from docs/source/processed_variable.rst
rename to docs/source/solvers/processed_variable.rst
index 9f1d62f2a7..4b4296061b 100644
--- a/docs/source/processed_variable.rst
+++ b/docs/source/solvers/processed_variable.rst
@@ -3,3 +3,6 @@ Post-Process Variables
 
 .. autoclass:: pybamm.ProcessedVariable
   :members:
+
+.. autoclass:: pybamm.ProcessedSymbolicVariable
+  :members:
diff --git a/docs/source/solvers/solution.rst b/docs/source/solvers/solution.rst
index 25a9977ec9..b39e1cae39 100644
--- a/docs/source/solvers/solution.rst
+++ b/docs/source/solvers/solution.rst
@@ -1,5 +1,5 @@
-Solution
-========
+Solutions
+=========
 
 .. autoclass:: pybamm._BaseSolution
   :members:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index aa25e2258d..121336d587 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -201,6 +201,8 @@ def version(formatted=False):
 # Solver classes
 #
 from .solvers.solution import Solution, _BaseSolution
+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
@@ -220,7 +222,6 @@ def version(formatted=False):
 #
 # other
 #
-from .processed_variable import ProcessedVariable
 from .quick_plot import QuickPlot, dynamic_plot, ax_min, ax_max
 
 from .simulation import Simulation, load_sim, is_notebook
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 6e0c480a33..5549b121d8 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -950,6 +950,14 @@ def _process_symbol(self, symbol):
 
             return new_symbol
 
+        elif isinstance(symbol, pybamm.InputParameter):
+            # Return a new copy of the input parameter, but set the expected size
+            # according to the domain of the input parameter
+            expected_size = self._get_variable_size(symbol)
+            new_input_parameter = symbol.new_copy()
+            new_input_parameter.set_expected_size(expected_size)
+            return new_input_parameter
+
         else:
             # Backup option: return new copy of the object
             try:
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index 8405fb76cb..e12b12565d 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -1,6 +1,7 @@
 #
 # Parameter classes
 #
+import numbers
 import numpy as np
 import pybamm
 
@@ -15,22 +16,32 @@ class InputParameter(pybamm.Symbol):
     ----------
     name : str
         name of the node
-
+    domain : iterable of str, or str
+        list of domains over which the node is valid (empty list indicates the symbol
+        is valid over all domains)
     """
 
-    def __init__(self, name):
-        super().__init__(name)
+    def __init__(self, name, domain=None):
+        # Expected shape defaults to 1
+        self._expected_size = 1
+        super().__init__(name, domain=domain)
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return InputParameter(self.name)
+        new_input_parameter = InputParameter(self.name, self.domain)
+        new_input_parameter._expected_size = self._expected_size
+        return new_input_parameter
+
+    def set_expected_size(self, size):
+        "Specify the size that the input parameter should be"
+        self._expected_size = size
 
     def _evaluate_for_shape(self):
         """
         Returns the scalar 'NaN' to represent the shape of a parameter.
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
         """
-        return np.nan
+        return np.nan * np.ones_like(self._expected_size)
 
     def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
@@ -44,10 +55,26 @@ def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
         if not isinstance(inputs, dict):
             # if the special input "shape test" is passed, just return 1
             if inputs == "shape test":
-                return 1
+                return np.ones_like(self._expected_size)
             raise TypeError("inputs should be a dictionary")
         try:
-            return inputs[self.name]
+            input_eval = inputs[self.name]
         # raise more informative error if can't find name in dict
         except KeyError:
             raise KeyError("Input parameter '{}' not found".format(self.name))
+
+        if isinstance(input_eval, numbers.Number):
+            input_shape = 1
+        else:
+            input_shape = input_eval.shape[0]
+        if input_shape == self._expected_size:
+            return input_eval
+        else:
+            raise ValueError(
+                "Input parameter '{}' was given an object of size '{}'".format(
+                    self.name, input_shape
+                )
+                + " but was expecting an object of size '{}'.".format(
+                    self._expected_size
+                )
+            )
diff --git a/pybamm/expression_tree/operations/jacobian.py b/pybamm/expression_tree/operations/jacobian.py
index b231c0c784..d714556615 100644
--- a/pybamm/expression_tree/operations/jacobian.py
+++ b/pybamm/expression_tree/operations/jacobian.py
@@ -76,7 +76,9 @@ def _jac(self, symbol, variable):
             jac = symbol._function_jac(children_jacs)
 
         elif isinstance(symbol, pybamm.Concatenation):
-            children_jacs = [child.jac(variable) for child in symbol.cached_children]
+            children_jacs = [
+                self.jac(child, variable) for child in symbol.cached_children
+            ]
             jac = symbol._concatenation_jac(children_jacs)
 
         else:
diff --git a/pybamm/expression_tree/operations/unpack_symbols.py b/pybamm/expression_tree/operations/unpack_symbols.py
index 2deed0d3c6..e28ecb4c81 100644
--- a/pybamm/expression_tree/operations/unpack_symbols.py
+++ b/pybamm/expression_tree/operations/unpack_symbols.py
@@ -22,7 +22,7 @@ def __init__(self, classes_to_find, unpacked_symbols=None):
 
     def unpack_list_of_symbols(self, list_of_symbols):
         """
-        Unpack a list of symbols. See :meth:`EquationUnpacker.unpack()`
+        Unpack a list of symbols. See :meth:`SymbolUnpacker.unpack()`
 
         Parameters
         ----------
@@ -65,7 +65,7 @@ def unpack_symbol(self, symbol):
             return unpacked
 
     def _unpack(self, symbol):
-        """ See :meth:`EquationUnpacker.unpack()`. """
+        """ See :meth:`SymbolUnpacker.unpack()`. """
 
         children = symbol.children
 
@@ -86,4 +86,3 @@ def _unpack(self, symbol):
                 child_vars = self.unpack_symbol(child)
                 found_vars.update(child_vars)
             return found_vars
-
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index cabcc24f48..bc30e48cf5 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -660,7 +660,9 @@ def evaluates_to_number(self):
         """
         result = self.evaluate_ignoring_errors()
 
-        if isinstance(result, numbers.Number):
+        if isinstance(result, numbers.Number) or (
+            isinstance(result, np.ndarray) and result.shape == ()
+        ):
             return True
         else:
             return False
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 516b937c89..131c5acd35 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -115,6 +115,7 @@ def __init__(self, name="Unnamed model"):
         self._jacobian = None
         self._jacobian_algebraic = None
         self.external_variables = []
+        self._input_parameters = None
 
         # Default behaviour is to use the jacobian and simplify
         self.use_jacobian = True
@@ -320,6 +321,25 @@ def timescale(self, value):
         "Set the timescale"
         self._timescale = value
 
+    @property
+    def input_parameters(self):
+        "Returns all the input parameters in the model"
+        if self._input_parameters is None:
+            self._input_parameters = self._find_input_parameters()
+        return self._input_parameters
+
+    def _find_input_parameters(self):
+        "Find all the input parameters in the model"
+        unpacker = pybamm.SymbolUnpacker(pybamm.InputParameter)
+        all_input_parameters = unpacker.unpack_list_of_symbols(
+            list(self.rhs.values())
+            + list(self.algebraic.values())
+            + list(self.initial_conditions.values())
+            + list(self.variables.values())
+            + [event.expression for event in self.events]
+        )
+        return list(all_input_parameters.values())
+
     def __getitem__(self, key):
         return self.rhs[key]
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 7db1e95860..d76e66dca2 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -183,15 +183,30 @@ def set_up(self, model, inputs=None):
             # set up Jacobian object, for re-use of dict
             jacobian = pybamm.Jacobian()
         else:
+            # Create placeholder inputs for evaluating rhs and algebraic sizes
+            placeholder_inputs = {}
+            for k, v in inputs.items():
+                if isinstance(v, casadi.MX):
+                    placeholder_inputs[k] = np.zeros(v.shape[0])
+                else:
+                    placeholder_inputs[k] = v
             # Convert model attributes to casadi
             t_casadi = casadi.MX.sym("t")
             y_diff = casadi.MX.sym(
                 "y_diff",
-                len(model.concatenated_rhs.evaluate(0, model.y0, inputs=inputs)),
+                len(
+                    model.concatenated_rhs.evaluate(
+                        0, model.y0, inputs=placeholder_inputs
+                    )
+                ),
             )
             y_alg = casadi.MX.sym(
                 "y_alg",
-                len(model.concatenated_algebraic.evaluate(0, model.y0, inputs=inputs)),
+                len(
+                    model.concatenated_algebraic.evaluate(
+                        0, model.y0, inputs=placeholder_inputs
+                    )
+                ),
             )
             y_casadi = casadi.vertcat(y_diff, y_alg)
             p_casadi = {}
@@ -575,7 +590,8 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         if (np.diff(t_eval) < 0).any():
             raise pybamm.SolverError("t_eval must increase monotonically")
 
-        # Non-dimensionalise t_eval
+        # Set up external variables and inputs
+        ext_and_inputs = self._set_up_ext_and_inputs(model, external_variables, inputs)
 
         # Make sure t_eval is monotonic
         if (np.diff(t_eval) < 0).any():
@@ -584,11 +600,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         # Set up
         timer = pybamm.Timer()
 
-        # Set up external variables and inputs
-        external_variables = external_variables or {}
-        inputs = inputs or {}
-        ext_and_inputs = {**external_variables, **inputs}
-
         # Set up (if not done already)
         if model not in self.models_set_up:
             self.set_up(model, ext_and_inputs)
@@ -699,13 +710,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         solution.set_up_time = set_up_time
         solution.solve_time = timer.time()
 
-        # Add model and inputs to solution
-        solution.model = model
-        solution.inputs = inputs
-
-        # Identify the event that caused termination
-        termination = self.get_termination_reason(solution, model.events)
-
         # Add model and inputs to solution
         solution.model = model
         solution.inputs = ext_and_inputs
@@ -880,6 +884,27 @@ def get_termination_reason(self, solution, events):
             solution.termination = "event: {}".format(termination_event)
             return "the termination event '{}' occurred".format(termination_event)
 
+    def _set_up_ext_and_inputs(self, model, external_variables, inputs):
+        "Set up external variables and input parameters"
+        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
+        for input_param in model.input_parameters:
+            name = input_param.name
+            if name not in inputs:
+                # Only allow symbolic inputs for CasadiAlgebraicSolver
+                if not isinstance(self, pybamm.CasadiAlgebraicSolver):
+                    raise pybamm.SolverError(
+                        "Only CasadiAlgebraicSolver can have symbolic inputs"
+                    )
+                inputs[name] = casadi.MX.sym(name, input_param._expected_size)
+
+        external_variables = external_variables or {}
+        ext_and_inputs = {**external_variables, **inputs}
+        return ext_and_inputs
+
 
 class SolverCallable:
     "A class that will be called by the solver when integrating"
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
index 8506f3a15e..70a5eda164 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -46,13 +46,20 @@ def _integrate(self, model, t_eval, inputs=None):
         t_eval : :class:`numpy.array`, size (k,)
             The times at which to compute the solution
         inputs : dict, optional
-            Any input parameters to pass to the model when solving
+            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.
         """
         y0 = model.y0
 
-        y = np.empty((len(y0), len(t_eval)))
+        y = None
 
         # Set up
+        # Record whether there are any symbolic inputs
+        has_symbolic_inputs = any(isinstance(v, casadi.MX) for v in inputs.values())
+
+        # Create casadi objects for the root-finder
         inputs = casadi.vertcat(*[x for x in inputs.values()])
         t_sym = casadi.MX.sym("t")
         y_sym = casadi.MX.sym("y_alg", y0.shape[0])
@@ -71,17 +78,23 @@ def _integrate(self, model, t_eval, inputs=None):
         for idx, t in enumerate(t_eval):
             # Evaluate algebraic with new t and previous y0, if it's already close
             # enough then keep it
-            if np.all(abs(model.algebraic_eval(t, y0, inputs)) < self.tol):
+            # We can't do this if there are symbolic inputs
+            if has_symbolic_inputs is False and np.all(
+                abs(model.casadi_algebraic(t, y0, inputs).full()) < self.tol
+            ):
                 pybamm.logger.debug(
                     "Keeping same solution at t={}".format(t * model.timescale_eval)
                 )
-                y[:, idx] = y0
-            # Otherwise calculate new y0
+                if y is None:
+                    y = y0
+                else:
+                    y = casadi.horzcat(y, y0)
+            # Otherwise calculate new y_sol
             else:
                 t_inputs = casadi.vertcat(t, inputs)
                 # Solve
                 try:
-                    y_sol = roots(y0, t_inputs).full().flatten()
+                    y_sol = roots(y0, t_inputs)
                     success = True
                     message = None
                     # Check final output
@@ -91,11 +104,18 @@ def _integrate(self, model, t_eval, inputs=None):
                     message = err.args[0]
                     fun = None
 
-                if success and np.all(casadi.fabs(fun) < self.tol):
+                # 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 np.all(casadi.fabs(fun) < self.tol)
+                ):
                     # update initial guess for the next iteration
                     y0 = y_sol
                     # update solution array
-                    y[:, idx] = y_sol
+                    if y is None:
+                        y = y_sol
+                    else:
+                        y = casadi.horzcat(y, y_sol)
                 elif not success:
                     raise pybamm.SolverError(
                         "Could not find acceptable solution: {}".format(message)
diff --git a/pybamm/solvers/processed_symbolic_variable.py b/pybamm/solvers/processed_symbolic_variable.py
new file mode 100644
index 0000000000..bd6e1f693f
--- /dev/null
+++ b/pybamm/solvers/processed_symbolic_variable.py
@@ -0,0 +1,230 @@
+#
+# 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.y.shape[0])
+        # Make all inputs symbolic first for converting to casadi
+        all_inputs_as_MX_dict = {}
+        symbolic_inputs_dict = {}
+        for key, value in solution.inputs.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()])
+        all_inputs = casadi.vertcat(*[p for p in solution.inputs.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_sol = solution.t
+        self.u_sol = solution.y
+        self.mesh = base_variable.mesh
+        self.symbolic_inputs_dict = symbolic_inputs_dict
+        self.symbolic_inputs_total_shape = symbolic_inputs.shape[0]
+        self.inputs = all_inputs
+        self.domain = base_variable.domain
+
+        self.base_eval = self.base_variable(solution.t[0], solution.y[:, 0], all_inputs)
+
+        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[0].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
+        for idx in range(len(self.t_sol)):
+            t = self.t_sol[idx]
+            u = self.u_sol[:, idx]
+            next_entries = self.base_variable(t, u, self.inputs)
+            if idx == 0:
+                entries = next_entries
+            else:
+                entries = casadi.horzcat(entries, next_entries)
+
+        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_sol)))
+
+        # Evaluate the base_variable index-by-index
+        for idx in range(len(self.t_sol)):
+            t = self.t_sol[idx]
+            u = self.u_sol[:, idx]
+            next_entries = self.base_variable(t, u, self.inputs)
+            if idx == 0:
+                entries = next_entries
+            else:
+                entries = casadi.horzcat(entries, next_entries)
+
+        # Get node values
+        nodes = self.mesh[0].nodes
+
+        # assign attributes for reference (either x_sol or r_sol)
+        self.entries = entries
+        self.dimensions = 1
+        if self.domain[0] in ["negative particle", "positive particle"]:
+            self.first_dimension = "r"
+            self.r_sol = nodes
+        elif self.domain[0] in [
+            "negative electrode",
+            "separator",
+            "positive electrode",
+        ]:
+            self.first_dimension = "x"
+            self.x_sol = nodes
+        elif self.domain == ["current collector"]:
+            self.first_dimension = "z"
+            self.z_sol = nodes
+        else:
+            self.first_dimension = "x"
+            self.x_sol = nodes
+
+        self.first_dim_pts = nodes
+
+    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.
+        """
+        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.
+        """
+        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))
+        # Check keys are consistent
+        if list(inputs_dict.keys()) != list(self.symbolic_inputs_dict.keys()):
+            raise ValueError(
+                "Inconsistent input keys: expected {}, actual {}".format(
+                    list(self.symbolic_inputs_dict.keys()), list(inputs_dict.keys())
+                )
+            )
+        inputs = casadi.vertcat(*[p for p in inputs_dict.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.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/processed_variable.py b/pybamm/solvers/processed_variable.py
similarity index 99%
rename from pybamm/processed_variable.py
rename to pybamm/solvers/processed_variable.py
index f10afd5042..e36ff46701 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/solvers/processed_variable.py
@@ -138,7 +138,7 @@ def initialise_1D(self):
             else:
                 entries[:, idx] = self.base_variable.evaluate(t, u, inputs=inputs)[:, 0]
 
-        # Process the discretisation to get x values
+        # Get node and edge values
         nodes = self.mesh[0].nodes
         edges = self.mesh[0].edges
         if entries.shape[0] == len(nodes):
diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py
index 725af3984c..097a25431d 100644
--- a/pybamm/solvers/solution.py
+++ b/pybamm/solvers/solution.py
@@ -1,6 +1,7 @@
 #
 # Solution class
 #
+import casadi
 import copy
 import numbers
 import numpy as np
@@ -40,6 +41,8 @@ def __init__(
         self, t, y, t_event=None, y_event=None, termination="final time", copy_this=None
     ):
         self._t = t
+        if isinstance(y, casadi.DM):
+            y = y.full()
         self._y = y
         self._t_event = t_event
         self._y_event = y_event
@@ -50,11 +53,13 @@ def __init__(
             self._model = None
             self.set_up_time = None
             self.solve_time = None
+            self.has_symbolic_inputs = False
         else:
             self._inputs = copy.copy(copy_this.inputs)
             self._model = copy_this.model
             self.set_up_time = copy_this.set_up_time
             self.solve_time = copy_this.solve_time
+            self.has_symbolic_inputs = copy_this.has_symbolic_inputs
 
         # initiaize empty variables and data
         self._variables = pybamm.FuzzyDict()
@@ -94,11 +99,18 @@ def inputs(self):
     @inputs.setter
     def inputs(self, inputs):
         "Updates the input values"
-        self._inputs = {}
-        for name, inp in inputs.items():
-            if isinstance(inp, numbers.Number):
-                inp = inp * np.ones_like(self.t)
-            self._inputs[name] = inp
+        # If there are symbolic inputs, just store them as given
+        if any(isinstance(v, casadi.MX) for v in inputs.values()):
+            self.has_symbolic_inputs = True
+            self._inputs = inputs
+        # Otherwise, make them the same size as the time vector
+        else:
+            self.has_symbolic_inputs = False
+            self._inputs = {}
+            for name, inp in inputs.items():
+                if isinstance(inp, numbers.Number):
+                    inp = inp * np.ones_like(self.t)
+                self._inputs[name] = inp
 
     @property
     def t_event(self):
@@ -142,13 +154,20 @@ def update(self, variables):
         # Process
         for key in variables:
             pybamm.logger.debug("Post-processing {}".format(key))
-            var = pybamm.ProcessedVariable(
-                self.model.variables[key], self, self._known_evals
-            )
+            # If there are symbolic inputs then we need to make a
+            # ProcessedSymbolicVariable
+            if self.has_symbolic_inputs is True:
+                var = pybamm.ProcessedSymbolicVariable(self.model.variables[key], self)
+
+            # Otherwise a standard ProcessedVariable is ok
+            else:
+                var = pybamm.ProcessedVariable(
+                    self.model.variables[key], self, self._known_evals
+                )
 
-            # Update known_evals in order to process any other variables faster
-            for t in var.known_evals:
-                self._known_evals[t].update(var.known_evals[t])
+                # Update known_evals in order to process any other variables faster
+                for t in var.known_evals:
+                    self._known_evals[t].update(var.known_evals[t])
 
             # Save variable and data
             self._variables[key] = var
@@ -248,6 +267,7 @@ class Solution(_BaseSolution):
 
     def __init__(self, t, y, t_event=None, y_event=None, termination="final time"):
         super().__init__(t, y, t_event, y_event, termination)
+        self.base_solution_class = _BaseSolution
 
     @property
     def sub_solutions(self):
@@ -278,7 +298,7 @@ def append(self, solution, start_index=1, create_sub_solutions=False):
         # functionality compared to normal solutions (can't append other solutions)
         if create_sub_solutions and not hasattr(self, "_sub_solutions"):
             self._sub_solutions = [
-                _BaseSolution(
+                self.base_solution_class(
                     self.t,
                     self.y,
                     self.t_event,
@@ -314,7 +334,7 @@ def append(self, solution, start_index=1, create_sub_solutions=False):
         # Append sub_solutions
         if create_sub_solutions:
             self._sub_solutions.append(
-                _BaseSolution(
+                self.base_solution_class(
                     solution.t,
                     solution.y,
                     solution.t_event,
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 5eb087d748..537d0d9a66 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -1189,6 +1189,18 @@ def test_mass_matirx_inverse(self):
             model.mass_matrix_inv.entries.toarray(), mass_inv.toarray()
         )
 
+    def test_process_input_variable(self):
+        disc = get_discretisation_for_testing()
+
+        a = pybamm.InputParameter("a")
+        a_disc = disc.process_symbol(a)
+        self.assertEqual(a_disc._expected_size, 1)
+
+        a = pybamm.InputParameter("a", ["negative electrode", "separator"])
+        a_disc = disc.process_symbol(a)
+        n = disc.mesh.combine_submeshes(*a.domain)[0].npts
+        self.assertEqual(a_disc._expected_size, n)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_expression_tree/test_input_parameter.py b/tests/unit/test_expression_tree/test_input_parameter.py
index 0df4ade271..f2d1e3b176 100644
--- a/tests/unit/test_expression_tree/test_input_parameter.py
+++ b/tests/unit/test_expression_tree/test_input_parameter.py
@@ -1,7 +1,7 @@
 #
 # Tests for the InputParameter class
 #
-import numbers
+import numpy as np
 import pybamm
 import unittest
 
@@ -13,9 +13,22 @@ def test_input_parameter_init(self):
         self.assertEqual(a.evaluate(inputs={"a": 1}), 1)
         self.assertEqual(a.evaluate(inputs={"a": 5}), 5)
 
+    def test_set_expected_size(self):
+        a = pybamm.InputParameter("a")
+        a.set_expected_size(10)
+        self.assertEqual(a._expected_size, 10)
+        y = np.linspace(0, 1, 10)
+        np.testing.assert_array_equal(a.evaluate(inputs={"a": y}), y)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Input parameter 'a' was given an object of size '1' but was expecting an "
+            "object of size '10'",
+        ):
+            a.evaluate(inputs={"a": 5})
+
     def test_evaluate_for_shape(self):
         a = pybamm.InputParameter("a")
-        self.assertIsInstance(a.evaluate_for_shape(), numbers.Number)
+        self.assertTrue(np.isnan(a.evaluate_for_shape()))
 
     def test_errors(self):
         a = pybamm.InputParameter("a")
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index 776ceebc05..1e19c0dfa9 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -107,6 +107,32 @@ def test_model_dict_behaviour(self):
         self.assertEqual(model[key], rhs[key])
         self.assertEqual(model[key], model.rhs[key])
 
+    def test_read_input_parameters(self):
+        # Read input parameters from different parts of the model
+        model = pybamm.BaseModel()
+        a = pybamm.InputParameter("a")
+        b = pybamm.InputParameter("b")
+        c = pybamm.InputParameter("c")
+        d = pybamm.InputParameter("d")
+        e = pybamm.InputParameter("e")
+        f = pybamm.InputParameter("f")
+
+        u = pybamm.Variable("u")
+        v = pybamm.Variable("v")
+        model.rhs = {u: -u * a}
+        model.algebraic = {v: v - b}
+        model.initial_conditions = {u: c, v: d}
+        model.events = [pybamm.Event("u=e", u - e)]
+        model.variables = {"v+f": v + f}
+
+        self.assertEqual(
+            set([x.name for x in model.input_parameters]),
+            set([x.name for x in [a, b, c, d, e, f]]),
+        )
+        self.assertTrue(
+            all(isinstance(x, pybamm.InputParameter) for x in model.input_parameters),
+        )
+
     def test_update(self):
         # model
         whole_cell = ["negative electrode", "separator", "positive electrode"]
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 3a7c2f130c..0540a874d2 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -50,6 +50,17 @@ def test_nonmonotonic_teval(self):
         ):
             solver.solve(model, np.array([1, 2, 3, 2]))
 
+    def test_block_symbolic_inputs(self):
+        solver = pybamm.BaseSolver(rtol=1e-2, atol=1e-4)
+        model = pybamm.BaseModel()
+        a = pybamm.Scalar(0)
+        p = pybamm.InputParameter("p")
+        model.rhs = {a: a * p}
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Only CasadiAlgebraicSolver can have symbolic inputs"
+        ):
+            solver.solve(model, np.array([1, 2, 3]))
+
     def test_ode_solver_fail_with_dae(self):
         model = pybamm.BaseModel()
         a = pybamm.Scalar(1)
diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index d69ccf92bd..efe1c0bab8 100644
--- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -5,6 +5,8 @@
 import pybamm
 import unittest
 import numpy as np
+from scipy.optimize import least_squares
+import tests
 
 
 class TestCasadiAlgebraicSolver(unittest.TestCase):
@@ -31,6 +33,23 @@ def test_simple_root_find(self):
         solution = solver.solve(model, np.linspace(0, 1, 10))
         np.testing.assert_array_equal(solution.y, -2)
 
+    def test_simple_root_find_correct_initial_guess(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + 2}
+        # initial guess gives right answer
+        model.initial_conditions = {var: -2}
+
+        # create discretisation
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, np.linspace(0, 1, 10))
+        np.testing.assert_array_equal(solution.y, -2)
+
     def test_root_find_fail(self):
         class Model:
             y0 = np.array([2])
@@ -88,7 +107,7 @@ def test_solve_with_input(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var")
         model = pybamm.BaseModel()
-        model.algebraic = {var: var + pybamm.InputParameter("value")}
+        model.algebraic = {var: var + pybamm.InputParameter("param")}
         model.initial_conditions = {var: 2}
 
         # create discretisation
@@ -97,10 +116,116 @@ def test_solve_with_input(self):
 
         # Solve
         solver = pybamm.CasadiAlgebraicSolver()
-        solution = solver.solve(model, np.linspace(0, 1, 10), inputs={"value": 7})
+        solution = solver.solve(model, np.linspace(0, 1, 10), inputs={"param": 7})
         np.testing.assert_array_equal(solution.y, -7)
 
 
+class TestCasadiAlgebraicSolverSensitivity(unittest.TestCase):
+    def test_solve_with_symbolic_input(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + pybamm.InputParameter("param")}
+        model.initial_conditions = {var: 2}
+        model.variables = {"var": var}
+
+        # create discretisation
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # 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)
+
+    def test_least_squares_fit(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var", domain="negative electrode")
+        model = pybamm.BaseModel()
+        p = pybamm.InputParameter("p")
+        q = pybamm.InputParameter("q")
+        model.algebraic = {var: (var - p)}
+        model.initial_conditions = {var: 3}
+        model.variables = {"objective": (var - q) ** 2 + (p - 3) ** 2}
+
+        # create discretisation
+        disc = tests.get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # 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()
+
+        # 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]})
+
+        # with jacobian
+        lsq_sol = least_squares(objective, [2, 2], jac=jac, method="lm")
+        np.testing.assert_array_almost_equal(lsq_sol.x, [3, 3], decimal=3)
+
+    def test_solve_with_symbolic_input_1D_scalar_input(self):
+        var = pybamm.Variable("var", "negative electrode")
+        model = pybamm.BaseModel()
+        param = pybamm.InputParameter("param")
+        model.algebraic = {var: var + param}
+        model.initial_conditions = {var: 2}
+        model.variables = {"var": var}
+
+        # create discretisation
+        disc = tests.get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # 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)
+
+    def test_solve_with_symbolic_input_1D_vector_input(self):
+        var = pybamm.Variable("var", "negative electrode")
+        model = pybamm.BaseModel()
+        param = pybamm.InputParameter("param", "negative electrode")
+        model.algebraic = {var: var + param}
+        model.initial_conditions = {var: 2}
+        model.variables = {"var": var}
+
+        # create discretisation
+        disc = tests.get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # Solve - scalar input
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, [0])
+        n = disc.mesh["negative electrode"][0].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
+        )
+        np.testing.assert_array_almost_equal(
+            solution["var"].value({"param": 2 * p}), -2 * p
+        )
+        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"].sensitivity({"param": p}), -np.eye(40)
+        )
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_solvers/test_processed_symbolic_variable.py b/tests/unit/test_solvers/test_processed_symbolic_variable.py
new file mode 100644
index 0000000000..cc4ad663f4
--- /dev/null
+++ b/tests/unit/test_solvers/test_processed_symbolic_variable.py
@@ -0,0 +1,317 @@
+#
+# 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)
+        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)
+        solution.inputs = {"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(ValueError, "Inconsistent input keys"):
+            processed_var.value({"not p": 3})
+        with self.assertRaisesRegex(ValueError, "Inconsistent input keys"):
+            processed_var.value({"q": 3, "p": 2})
+
+    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)
+        solution.inputs = {"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 = [0]
+        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        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)
+        processed_eqn = pybamm.ProcessedSymbolicVariable(eqn_sol, sol)
+        np.testing.assert_array_equal(
+            processed_eqn.value(), y_sol + x_sol[:, np.newaxis]
+        )
+
+    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 = [0]
+        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
+
+        sol = pybamm.Solution(t_sol, y_sol)
+        sol.inputs = {"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)
+        sol.inputs = {"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": -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"])
+        p.set_expected_size(65)
+        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 = [0]
+        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        sol.inputs = {"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 = [0]
+        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        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 = [0]
+        y_sol = np.ones_like(r_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        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 = [0]
+        y_sol = np.ones_like(z_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        pybamm.ProcessedSymbolicVariable(var_sol, sol)
+
+        # Other domain
+        var = pybamm.Variable("var", domain=["line"])
+        x = pybamm.SpatialVariable("x", domain=["line"])
+
+        geometry = pybamm.Geometry(
+            {
+                "line": {
+                    "primary": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
+                }
+            }
+        )
+        submesh_types = {"line": pybamm.MeshGenerator(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 = [0]
+        y_sol = np.ones_like(x_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        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 = [0]
+        y_sol = np.ones_like(r_sol)[:, np.newaxis] * 5
+        sol = pybamm.Solution(t_sol, y_sol)
+        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_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py
similarity index 99%
rename from tests/unit/test_processed_variable.py
rename to tests/unit/test_solvers/test_processed_variable.py
index 09e102987b..224fe24d22 100644
--- a/tests/unit/test_processed_variable.py
+++ b/tests/unit/test_solvers/test_processed_variable.py
@@ -32,7 +32,6 @@ def test_processed_variable_1D(self):
         x_sol = disc.process_symbol(x).entries[:, 0]
         var_sol = disc.process_symbol(var)
         eqn_sol = disc.process_symbol(eqn)
-        eqn_sol.mesh = disc.mesh.combine_submeshes(*eqn.domain)
         t_sol = np.linspace(0, 1)
         y_sol = np.ones_like(x_sol)[:, np.newaxis] * np.linspace(0, 5)
 
@@ -223,7 +222,6 @@ def test_processed_var_1D_interpolation(self):
         x_sol = disc.process_symbol(x).entries[:, 0]
         var_sol = disc.process_symbol(var)
         eqn_sol = disc.process_symbol(eqn)
-        eqn_sol.mesh = disc.mesh.combine_submeshes(*eqn.domain)
         t_sol = np.linspace(0, 1)
         y_sol = x_sol[:, np.newaxis] * np.linspace(0, 5)