diff --git a/CHANGELOG.md b/CHANGELOG.md
index ca6f825ebb..b3a3a23a06 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added `InputParameter` node for quickly changing parameter values ([#752](https://github.com/pybamm-team/PyBaMM/pull/752))
 -   Generalized importing of external variables ([#728](https://github.com/pybamm-team/PyBaMM/pull/728))
 -   Separated active and inactive material volume fractions ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
 -   Added submodels for tortuosity ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
diff --git a/docs/source/expression_tree/index.rst b/docs/source/expression_tree/index.rst
index 3eedd7b622..3efd97fb3f 100644
--- a/docs/source/expression_tree/index.rst
+++ b/docs/source/expression_tree/index.rst
@@ -16,5 +16,6 @@ Expression Tree
   concatenations
   broadcasts
   functions
+  input_parameter
   interpolant
   operations/index
diff --git a/docs/source/expression_tree/input_parameter.rst b/docs/source/expression_tree/input_parameter.rst
new file mode 100644
index 0000000000..ebbc13e7e5
--- /dev/null
+++ b/docs/source/expression_tree/input_parameter.rst
@@ -0,0 +1,5 @@
+Input Parameter
+===============
+
+.. autoclass:: pybamm.InputParameter
+  :members:
diff --git a/docs/tutorials/add-parameter-values.rst b/docs/tutorials/add-parameter-values.rst
index 3cc8ed090e..9472e9a872 100644
--- a/docs/tutorials/add-parameter-values.rst
+++ b/docs/tutorials/add-parameter-values.rst
@@ -21,56 +21,36 @@ For an example of how the parameter values work, see the
 Adding a set of parameters values
 ---------------------------------
 
-There are two ways to add parameter sets:
+Parameter sets are split by material into ``anodes``, ``separators``, ``cathodes``, ``electrolytes``, ``cells`` (for cell geometries and thermal properties) and ``experiments`` (for initial conditions and charge/discharge rates).
+To add a new parameter set in one of these subcategories, first create a new folder in the appropriate chemistry folder: for example, to add a new anode chemistry for lithium-ion, add a subfolder ``input/parameters/lithium-ion/anodes/new_anode_chemistry_AuthorYear``. 
+This subfolder should then contain:
 
-1. **Complete parameter file**: Parameter sets should be added as csv files in the appropriate chemistry folder in ``input/parameters/`` (add a new folder if no parameters exist for that chemistry yet).
-The expected structure of the csv file is
+- a csv file ``parameters.csv`` with all the relevant scalar parameters. The expected structure of the csv file is:
 
 +-------------------------+----------------------+-----------------------+-------------+
 | Name [Units]            | Value                | Reference             | Notes       |
 +=========================+======================+=======================+=============+
-| Example [m]             | 13                   | bloggs2019            | an example  |
+| Example [m]             | 13                   | AuthorYear            | an example  |
 +-------------------------+----------------------+-----------------------+-------------+
 
 Empty lines, and lines starting with ``#``, will be ignored.
 
-2. **Parameters for a single material**: It's possible to add parameters for a single material (e.g. anode) and then re-use existing parameters for the other materials. To do this, add a new subfolder with a ``README.md`` for references and csv file of parameters (e.g. ``input/parameters/lithium-ion/anodes/new_anode_chemistry_Bloggs2019/``). Then this file can be referenced using the ``chemistry`` option in ``ParameterValues``, e.g.
+- a ``README.md`` file with information on where these parameters came from
+- python files for any functions, which should be referenced from the ``parameters.csv`` file (see ``Adding a Function`` below)
+- csv files for any data to be interpolated, which should be referenced from the ``parameters.csv`` file (see ``Adding data for interpolation`` below)
 
-.. code-block:: python
-
-    param = pybamm.ParameterValues(
-        chemistry={
-            "chemistry": "lithium-ion",
-            "cell": "kokam_Marquis2019",
-            "anode": "new_anode_chemistry_Bloggs2019",
-            "separator": "separator_Marquis2019",
-            "cathode": "lico2_Marquis2019",
-            "electrolyte": "lipf6_Marquis2019",
-            "experiment": "1C_discharge_from_full_Marquis2019",
-        }
-    )
-
-or, equivalently in this case (since the only difference from the standard parameters from Marquis et al. is the set of anode parameters),
-
-.. code-block:: python
-
-    param = pybamm.ParameterValues(
-        chemistry={
-            **pybamm.parameter_sets.Marquis2019,
-            "anode": "new_anode_chemistry_Bloggs2019",
-        }
-    )
+The easiest way to start is to copy an existing file (e.g. ````input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019``) and replace all entries in all files as appropriate
 
 Adding a function
 -----------------
 
 Functions should be added as Python functions under a file with the same name in the appropriate chemistry folder in ``input/parameters/``.
 These Python functions should be documented with references explaining where they were obtained.
-For example, we would put the following Python function in a file ``input/parameters/lead-acid/diffusivity_Bloggs2019.py``
+For example, we would put the following Python function in a file ``input/parameters/lithium_ion/anodes/new_anode_chemistry_AuthorYear/diffusivity_AuthorYear.py``
 
 .. code-block:: python
 
-    def diffusivity_Bloggs2019(c_e):
+    def diffusivity_AuthorYear(c_e):
         """
         Dimensional Fickian diffusivity in the electrolyte [m2.s-1], from [1]_, as a
         function of the electrolyte concentration c_e [mol.m-3].
@@ -89,10 +69,92 @@ called (must be in the same folder), with the tag ``[function]``, for example:
 +---------------------+--------------------------------------+--------------+-------------+
 | Name [Units]        | Value                                |  Reference   | Notes       |
 +=====================+======================================+==============+=============+
-| Example [m2.s-1]    | [function]diffusivity_Bloggs2019     | bloggs2019   | a function  |
+| Example [m2.s-1]    | [function]diffusivity_AuthorYear     | AuthorYear   | a function  |
 +---------------------+--------------------------------------+--------------+-------------+
 
+Adding data for interpolation
+-----------------------------
+
+Data should be added as as csv file in the appropriate chemistry folder in ``input/parameters/``.
+For example, we would put the following data in a file ``input/parameters/lithium_ion/anodes/new_anode_chemistry_AuthorYear/diffusivity_AuthorYear.csv``
+
++--------------------------+--------------------------+
+| # concentration [mol/m3] | Diffusivity [m2/s]       |
++==========================+==========================+
+| 0.000000000000000000e+00 | 4.714135898019971016e+00 |
+| 2.040816326530612082e-02 | 4.708899441575220557e+00 |
+| 4.081632653061224164e-02 | 4.702448345762175741e+00 |
+| 6.122448979591836593e-02 | 4.694558534379876136e+00 |
+| 8.163265306122448328e-02 | 4.684994372928071193e+00 |
+| 1.020408163265306006e-01 | 4.673523893805322516e+00 |
+| 1.224489795918367319e-01 | 4.659941254449398329e+00 |
+| 1.428571428571428492e-01 | 4.644096031712390271e+00 |
++--------------------------+--------------------------+
+
+Empty lines, and lines starting with ``#``, will be ignored.
+
+Then, this data should be added to the parameter file from which it will be
+called (must be in the same folder), with the tag ``[data]``, for example:
+
++---------------------+----------------------------------+--------------+-------------+
+| Name [Units]        | Value                            |  Reference   | Notes       |
++=====================+==================================+==============+=============+
+| Example [m2.s-1]    | [data]diffusivity_AuthorYear     | AuthorYear   | some data   |
++---------------------+----------------------------------+--------------+-------------+
+
+Using new parameters
+--------------------
+
+If you have added a whole new set of parameters, then you can create a new parameter set in ``pybamm/parameters/parameter_sets.py``, by just adding a new dictionary to that file, for example
+
+.. code-block:: python
+
+    AuthorYear = {
+        "chemistry": "lithium-ion",
+        "cell": "new_cell_AuthorYear",
+        "anode": "new_anode_AuthorYear",
+        "separator": "new_separator_AuthorYear",
+        "cathode": "new_cathode_AuthorYear",
+        "electrolyte": "new_electrolyte_AuthorYear",
+        "experiment": "new_experiment_AuthorYear",
+    }
+
+Then, to use these new parameters, use:
+
+.. code-block:: python
+
+    param = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.AuthorYear)
+
+Note that you can re-use existing parameter subsets instead of creating new ones (for example, you could just replace "experiment": "new_experiment_AuthorYear" with "experiment": "1C_discharge_from_full_Marquis2019" in the above dictionary).
+
+It's also possible to add parameters for a single material (e.g. anode) and then re-use existing parameters for the other materials, without adding a parameter set to ``pybamm/parameters/parameter_sets.py``.
+
+.. code-block:: python
+
+    param = pybamm.ParameterValues(
+        chemistry={
+            "chemistry": "lithium-ion",
+            "cell": "kokam_Marquis2019",
+            "anode": "new_anode_chemistry_AuthorYear",
+            "separator": "separator_Marquis2019",
+            "cathode": "lico2_Marquis2019",
+            "electrolyte": "lipf6_Marquis2019",
+            "experiment": "1C_discharge_from_full_Marquis2019",
+        }
+    )
+
+or, equivalently in this case (since the only difference from the standard parameters from Marquis et al. is the set of anode parameters),
+
+.. code-block:: python
+
+    param = pybamm.ParameterValues(
+        chemistry={
+            **pybamm.parameter_sets.Marquis2019,
+            "anode": "new_anode_chemistry_AuthorYear",
+        }
+    )
 
+See the `"Getting Started" tutorial <https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/Getting%20Started/Tutorial%202%20-%20Setting%20Parameter%20Values.ipynb>`_ for examples of setting parameters in action.
 
 Unit tests for the new class
 ----------------------------
@@ -105,13 +167,13 @@ Test on the models
 
 In theory, any existing model can now be solved using the new parameters instead of their default parameters, with no extra work from here.
 To test this, add something like the following test to one of the model test files
-(e.g. `DFN <https://github.com/pybamm-team/PyBaMM/blob/master/tests/test_models/test_lithium_ion/test_lithium_ion_dfn.py>`_):
+(e.g. `DFN <https://github.com/pybamm-team/PyBaMM/blob/master/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py>`_):
 
 .. code-block:: python
 
     def test_my_new_parameters(self):
         model = pybamm.lithium_ion.DFN()
-        parameter_values = pybamm.ParameterValues("path/to/parameter/file.csv")
+        parameter_values = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.AuthorYear)
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all()
 
diff --git a/docs/tutorials/add-solver.rst b/docs/tutorials/add-solver.rst
index 1ee8e4f328..fbdf60bcd3 100644
--- a/docs/tutorials/add-solver.rst
+++ b/docs/tutorials/add-solver.rst
@@ -74,7 +74,7 @@ Test on the models
 
 In theory, any existing model can now be solved using `MyFastDaeSolver` instead of their default solvers, with no extra work from here.
 To test this, add something like the following test to one of the model test files
-(e.g. `DFN <https://github.com/pybamm-team/PyBaMM/blob/master/tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py>`_):
+(e.g. `DFN <https://github.com/pybamm-team/PyBaMM/blob/master/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py>`_):
 
 .. code-block:: python
 
diff --git a/docs/tutorials/add-spatial-method.rst b/docs/tutorials/add-spatial-method.rst
index 2d6ad71c31..00cf2ae928 100644
--- a/docs/tutorials/add-spatial-method.rst
+++ b/docs/tutorials/add-spatial-method.rst
@@ -87,7 +87,7 @@ Test on the models
 
 In theory, any existing model can now be discretised using ``MyFastMethod`` instead of their default spatial methods, with no extra work from here.
 To test this, add something like the following test to one of the model test files
-(e.g. `DFN <https://github.com/pybamm-team/PyBaMM/blob/master/tests/unit/test_models/test_lithium_ion/test_lithium_ion_dfn.py>`_):
+(e.g. `DFN <https://github.com/pybamm-team/PyBaMM/blob/master/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py>`_):
 
 .. code-block:: python
 
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 6db2ab6199..a97a49f57d 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -98,6 +98,7 @@ def version(formatted=False):
 from .expression_tree.unary_operators import *
 from .expression_tree.functions import *
 from .expression_tree.interpolant import Interpolant
+from .expression_tree.input_parameter import InputParameter
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.broadcasts import (
     Broadcast,
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index eafece3c57..dffbbaa367 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -174,21 +174,21 @@ def _binary_new_copy(self, left, right):
         "Default behaviour for new_copy"
         return self.__class__(left, right)
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             id = self.id
             try:
                 return known_evals[id], known_evals
             except KeyError:
-                left, known_evals = self.left.evaluate(t, y, known_evals)
-                right, known_evals = self.right.evaluate(t, y, known_evals)
+                left, known_evals = self.left.evaluate(t, y, u, known_evals)
+                right, known_evals = self.right.evaluate(t, y, u, known_evals)
                 value = self._binary_evaluate(left, right)
                 known_evals[id] = value
                 return value, known_evals
         else:
-            left = self.left.evaluate(t, y)
-            right = self.right.evaluate(t, y)
+            left = self.left.evaluate(t, y, u)
+            right = self.right.evaluate(t, y, u)
             return self._binary_evaluate(left, right)
 
     def evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index c90d2874e0..2508b4cc8b 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -52,20 +52,22 @@ def _concatenation_evaluate(self, children_eval):
         else:
             return self.concatenation_function(children_eval)
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         children = self.cached_children
         if known_evals is not None:
             if self.id not in known_evals:
                 children_eval = [None] * len(children)
                 for idx, child in enumerate(children):
-                    children_eval[idx], known_evals = child.evaluate(t, y, known_evals)
+                    children_eval[idx], known_evals = child.evaluate(
+                        t, y, u, known_evals
+                    )
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
             children_eval = [None] * len(children)
             for idx, child in enumerate(children):
-                children_eval[idx] = child.evaluate(t, y)
+                children_eval[idx] = child.evaluate(t, y, u)
             return self._concatenation_evaluate(children_eval)
 
     def new_copy(self):
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index 8ce840815f..0014b58fbb 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -14,8 +14,8 @@ class Function(pybamm.Symbol):
     ----------
     function : method
         A function can have 0 or many inputs. If no inputs are given, self.evaluate()
-        simply returns func(). Otherwise, self.evaluate(t, y) returns
-        func(child0.evaluate(t, y), child1.evaluate(t, y), etc).
+        simply returns func(). Otherwise, self.evaluate(t, y, u) returns
+        func(child0.evaluate(t, y, u), child1.evaluate(t, y, u), etc).
     children : :class:`pybamm.Symbol`
         The children nodes to apply the function to
     derivative : str, optional
@@ -152,19 +152,19 @@ def _function_jac(self, children_jacs):
 
         return jacobian
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
                 evaluated_children = [None] * len(self.children)
                 for i, child in enumerate(self.children):
                     evaluated_children[i], known_evals = child.evaluate(
-                        t, y, known_evals
+                        t, y, known_evals=known_evals
                     )
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
         else:
-            evaluated_children = [child.evaluate(t, y) for child in self.children]
+            evaluated_children = [child.evaluate(t, y, u) for child in self.children]
             return self._function_evaluate(evaluated_children)
 
     def evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
new file mode 100644
index 0000000000..67b8a87551
--- /dev/null
+++ b/pybamm/expression_tree/input_parameter.py
@@ -0,0 +1,53 @@
+#
+# Parameter classes
+#
+import numpy as np
+import pybamm
+
+
+class InputParameter(pybamm.Symbol):
+    """A node in the expression tree representing an input parameter
+
+    This node's value can be set at the point of solving, allowing parameter estimation
+    and control
+
+    Parameters
+    ----------
+    name : str
+        name of the node
+
+    """
+
+    def __init__(self, name):
+        super().__init__(name)
+
+    def new_copy(self):
+        """ See :meth:`pybamm.Symbol.new_copy()`. """
+        return InputParameter(self.name)
+
+    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
+
+    def _jac(self, variable):
+        """ See :meth:`pybamm.Symbol._jac()`. """
+        return pybamm.Scalar(0)
+
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+        # u should be a dictionary
+        # convert 'None' to empty dictionary for more informative error
+        if u is None:
+            u = {}
+        if not isinstance(u, dict):
+            # if the special input "shape test" is passed, just return 1
+            if u == "shape test":
+                return 1
+            raise TypeError("inputs u should be a dictionary")
+        try:
+            return u[self.name]
+        # raise more informative error if can't find name in dict
+        except KeyError:
+            raise KeyError("Input parameter '{}' not found".format(self.name))
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index e6e0758410..983008d706 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -11,36 +11,43 @@ class CasadiConverter(object):
     def __init__(self, casadi_symbols=None):
         self._casadi_symbols = casadi_symbols or {}
 
-    def convert(self, symbol, t=None, y=None):
+    def convert(self, symbol, t=None, y=None, u=None):
         """
-        This function recurses down the tree, applying any simplifications defined in
-        classes derived from pybamm.Symbol. E.g. any expression multiplied by a
-        pybamm.Scalar(0) will be simplified to a pybamm.Scalar(0).
-        If a symbol has already been simplified, the stored value is returned.
+        This function recurses down the tree, converting the PyBaMM expression tree to
+        a CasADi expression tree
 
         Parameters
         ----------
         symbol : :class:`pybamm.Symbol`
             The symbol to convert
+        t : :class:`casadi.MX`
+            A casadi symbol representing time
+        y : :class:`casadi.MX`
+            A casadi symbol representing state vectors
+        u : dict
+            A dictionary of casadi symbols representing inputs
 
         Returns
         -------
-        CasADi symbol
-            The convert symbol
+        :class:`casadi.MX`
+            The converted symbol
         """
-
         try:
             return self._casadi_symbols[symbol.id]
         except KeyError:
-            casadi_symbol = self._convert(symbol, t, y)
+            # Change u to empty dictionary if it's None
+            u = u or {}
+            casadi_symbol = self._convert(symbol, t, y, u)
             self._casadi_symbols[symbol.id] = casadi_symbol
 
             return casadi_symbol
 
-    def _convert(self, symbol, t, y):
+    def _convert(self, symbol, t=None, y=None, u=None):
         """ See :meth:`CasadiConverter.convert()`. """
-        if isinstance(symbol, (pybamm.Scalar, pybamm.Array, pybamm.Time)):
-            return casadi.MX(symbol.evaluate(t, y))
+        if isinstance(
+            symbol, (pybamm.Scalar, pybamm.Array, pybamm.Time, pybamm.InputParameter)
+        ):
+            return casadi.MX(symbol.evaluate(t, y, u))
 
         elif isinstance(symbol, pybamm.StateVector):
             if y is None:
@@ -50,8 +57,8 @@ def _convert(self, symbol, t, y):
         elif isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
             # process children
-            converted_left = self.convert(left, t, y)
-            converted_right = self.convert(right, t, y)
+            converted_left = self.convert(left, t, y, u)
+            converted_right = self.convert(right, t, y, u)
             if isinstance(symbol, pybamm.Outer):
                 return casadi.kron(converted_left, converted_right)
             else:
@@ -59,14 +66,14 @@ def _convert(self, symbol, t, y):
                 return symbol._binary_evaluate(converted_left, converted_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
-            converted_child = self.convert(symbol.child, t, y)
+            converted_child = self.convert(symbol.child, t, y, u)
             if isinstance(symbol, pybamm.AbsoluteValue):
                 return casadi.fabs(converted_child)
             return symbol._unary_evaluate(converted_child)
 
         elif isinstance(symbol, pybamm.Function):
             converted_children = [
-                self.convert(child, t, y) for child in symbol.children
+                self.convert(child, t, y, u) for child in symbol.children
             ]
             # Special functions
             if symbol.function == np.min:
@@ -97,7 +104,7 @@ def _convert(self, symbol, t, y):
                 return symbol._function_evaluate(converted_children)
         elif isinstance(symbol, pybamm.Concatenation):
             converted_children = [
-                self.convert(child, t, y) for child in symbol.children
+                self.convert(child, t, y, u) for child in symbol.children
             ]
             if isinstance(symbol, (pybamm.NumpyConcatenation, pybamm.SparseStack)):
                 return casadi.vertcat(*converted_children)
diff --git a/pybamm/expression_tree/operations/evaluate.py b/pybamm/expression_tree/operations/evaluate.py
index 43065a91b8..6558357c33 100644
--- a/pybamm/expression_tree/operations/evaluate.py
+++ b/pybamm/expression_tree/operations/evaluate.py
@@ -175,6 +175,9 @@ def find_symbols(symbol, constant_symbols, variable_symbols):
     elif isinstance(symbol, pybamm.Time):
         symbol_str = "t"
 
+    elif isinstance(symbol, pybamm.InputParameter):
+        symbol_str = "u['{}']".format(symbol.name)
+
     else:
         raise NotImplementedError(
             "Not implemented for a symbol of type '{}'".format(type(symbol))
@@ -270,7 +273,7 @@ def __init__(self, symbol):
             self._result_var, "return" + self._result_var, "eval"
         )
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """
         Acts as a drop-in replacement for :func:`pybamm.Symbol.evaluate`
         """
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index a5af685ab2..ca902f5172 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -459,7 +459,7 @@ def _base_evaluate(self, t=None, y=None):
             )
         )
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """Evaluate expression tree (wrapper to allow using dict of known values).
         If the dict 'known_evals' is provided, the dict is searched for self.id; if
         self.id is in the keys, return that value; otherwise, evaluate using
@@ -498,6 +498,7 @@ def evaluate_for_shape(self):
 
     def is_constant(self):
         """returns true if evaluating the expression is not dependent on `t` or `y`
+        or `u`
 
         See Also
         --------
@@ -505,8 +506,13 @@ def is_constant(self):
 
         """
         # if any of the nodes are instances of any of these types, then the whole
-        # expression depends on either t or y
-        search_types = (pybamm.Variable, pybamm.StateVector, pybamm.IndependentVariable)
+        # expression depends on either t or y or u
+        search_types = (
+            pybamm.Variable,
+            pybamm.StateVector,
+            pybamm.Time,
+            pybamm.InputParameter,
+        )
 
         # do the search, return true if no relevent nodes are found
         return not any((isinstance(n, search_types)) for n in self.pre_order())
@@ -514,8 +520,8 @@ def is_constant(self):
     def evaluate_ignoring_errors(self):
         """
         Evaluates the expression. If a node exists in the tree that cannot be evaluated
-        as a scalar or vectr (e.g. Parameter, Variable, StateVector), then None is
-        returned. Otherwise the result of the evaluation is given
+        as a scalar or vector (e.g. Parameter, Variable, StateVector, InputParameter),
+        then None is returned. Otherwise the result of the evaluation is given
 
         See Also
         --------
@@ -523,19 +529,18 @@ def evaluate_ignoring_errors(self):
 
         """
         try:
-            result = self.evaluate(t=0)
+            result = self.evaluate(t=0, u="shape test")
         except NotImplementedError:
-            # return false if NotImplementedError is raised
+            # return None if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
             return None
         except TypeError as error:
-            # return false if specific TypeError is raised
+            # return None if specific TypeError is raised
             # (there is a e.g. StateVector in the tree)
             if error.args[0] == "StateVector cannot evaluate input 'y=None'":
                 return None
             else:
                 raise error
-
         return result
 
     def evaluates_to_number(self):
@@ -579,12 +584,12 @@ def simplify(self, simplified_symbols=None):
         """ Simplify the expression tree. See :class:`pybamm.Simplification`. """
         return pybamm.Simplification(simplified_symbols).simplify(self)
 
-    def to_casadi(self, t=None, y=None, casadi_symbols=None):
+    def to_casadi(self, t=None, y=None, u=None, casadi_symbols=None):
         """
         Convert the expression tree to a CasADi expression tree.
         See :class:`pybamm.CasadiConverter`.
         """
-        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y)
+        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y, u)
 
     def new_copy(self):
         """
@@ -614,7 +619,7 @@ def shape(self):
         # Try with some large y, to avoid having to use pre_order (slow)
         try:
             y = np.linspace(0.1, 0.9, int(1e4))
-            evaluated_self = self.evaluate(0, y)
+            evaluated_self = self.evaluate(0, y, u="shape test")
         # If that fails, fall back to calculating how big y should really be
         except ValueError:
             state_vectors_in_node = [
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 5ce164498a..0fcd87c066 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -63,15 +63,15 @@ def _unary_evaluate(self, child):
         """Perform unary operation on a child. """
         raise NotImplementedError
 
-    def evaluate(self, t=None, y=None, known_evals=None):
+    def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
-                child, known_evals = self.child.evaluate(t, y, known_evals)
+                child, known_evals = self.child.evaluate(t, y, u, known_evals)
                 known_evals[self.id] = self._unary_evaluate(child)
             return known_evals[self.id], known_evals
         else:
-            child = self.child.evaluate(t, y)
+            child = self.child.evaluate(t, y, u)
             return self._unary_evaluate(child)
 
     def evaluate_for_shape(self):
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index b16e88a637..50b2bb35a6 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -156,7 +156,7 @@ def options(self, extra_options):
             "particle": "Fickian diffusion",
             "thermal": "isothermal",
             "thermal current collector": False,
-            "external submodels": []
+            "external submodels": [],
         }
         options = default_options
         # any extra options overwrite the default options
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 69dd1f6650..26c7b51b1d 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -175,6 +175,8 @@ def update(self, values, check_conflict=False, path=""):
                     # Special case (hacky) for zero current
                     elif value == "[zero]":
                         super().__setitem__(name, 0)
+                    elif value == "[input]":
+                        super().__setitem__(name, pybamm.InputParameter(name))
                     # Anything else should be a converted to a float
                     else:
                         super().__setitem__(name, float(value))
@@ -419,8 +421,13 @@ def _process_symbol(self, symbol):
 
         if isinstance(symbol, pybamm.Parameter):
             value = self[symbol.name]
-            # Scalar inherits name (for updating parameters) and domain (for Broadcast)
-            return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
+            if isinstance(value, numbers.Number):
+                # Scalar inherits name (for updating parameters) and domain (for
+                # Broadcast)
+                return pybamm.Scalar(value, name=symbol.name, domain=symbol.domain)
+            elif isinstance(value, pybamm.InputParameter):
+                value.domain = symbol.domain
+                return value
 
         elif isinstance(symbol, pybamm.FunctionParameter):
             new_children = [self.process_symbol(child) for child in symbol.children]
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index a4db4a6462..5081276be6 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -7,7 +7,9 @@
 import scipy.interpolate as interp
 
 
-def post_process_variables(variables, t_sol, u_sol, mesh=None, interp_kind="linear"):
+def post_process_variables(
+    variables, t_sol, u_sol, mesh=None, inputs=None, interp_kind="linear"
+):
     """
     Post-process all variables in a model
 
@@ -23,6 +25,8 @@ def post_process_variables(variables, t_sol, u_sol, mesh=None, interp_kind="line
     mesh : :class:`pybamm.Mesh`
         The mesh used to solve, used here to calculate the reference x values for
         interpolation
+    inputs : dict, optional
+        Any input parameters to pass to the model
     interp_kind : str
         The method to use for interpolation
 
@@ -36,7 +40,7 @@ def post_process_variables(variables, t_sol, u_sol, mesh=None, interp_kind="line
     for var, eqn in variables.items():
         pybamm.logger.debug("Post-processing {}".format(var))
         processed_variables[var] = ProcessedVariable(
-            eqn, t_sol, u_sol, mesh, interp_kind, known_evals
+            eqn, t_sol, u_sol, mesh, inputs, interp_kind, known_evals
         )
 
         for t in known_evals:
@@ -64,6 +68,8 @@ class ProcessedVariable(object):
     mesh : :class:`pybamm.Mesh`
         The mesh used to solve, used here to calculate the reference x values for
         interpolation
+    inputs : dict, optional
+            Any input parameters to pass to the model
     interp_kind : str
         The method to use for interpolation
     """
@@ -74,6 +80,7 @@ def __init__(
         t_sol,
         u_sol,
         mesh=None,
+        inputs=None,
         interp_kind="linear",
         known_evals=None,
     ):
@@ -81,6 +88,7 @@ def __init__(
         self.t_sol = t_sol
         self.u_sol = u_sol
         self.mesh = mesh
+        self.inputs = inputs or {}
         self.interp_kind = interp_kind
         self.domain = base_variable.domain
         self.auxiliary_domains = base_variable.auxiliary_domains
@@ -88,7 +96,10 @@ def __init__(
 
         if self.known_evals:
             self.base_eval, self.known_evals[t_sol[0]] = base_variable.evaluate(
-                t_sol[0], u_sol[:, 0], self.known_evals[t_sol[0]]
+                t_sol[0],
+                u_sol[:, 0],
+                self.inputs,
+                known_evals=self.known_evals[t_sol[0]],
             )
         else:
             self.base_eval = base_variable.evaluate(t_sol[0], u_sol[:, 0])
@@ -131,10 +142,12 @@ def initialise_1D(self):
             t = self.t_sol[idx]
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
-                    t, self.u_sol[:, idx], self.known_evals[t]
+                    t, self.u_sol[:, idx], self.inputs, known_evals=self.known_evals[t]
                 )
             else:
-                entries[idx] = self.base_variable.evaluate(t, self.u_sol[:, idx])
+                entries[idx] = self.base_variable.evaluate(
+                    t, self.u_sol[:, idx], self.inputs
+                )
 
         # No discretisation provided, or variable has no domain (function of t only)
         self._interpolation_function = interp.interp1d(
@@ -158,12 +171,12 @@ def initialise_2D(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, idx] = eval_and_known_evals[0][:, 0]
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
-                entries[:, idx] = self.base_variable.evaluate(t, u)[:, 0]
+                entries[:, idx] = self.base_variable.evaluate(t, u, self.inputs)[:, 0]
 
         # Process the discretisation to get x values
         nodes = self.mesh.combine_submeshes(*self.domain)[0].nodes
@@ -301,7 +314,7 @@ def initialise_3D(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(
                     eval_and_known_evals[0],
@@ -311,7 +324,7 @@ def initialise_3D(self):
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u),
+                    self.base_variable.evaluate(t, u, self.inputs),
                     [first_dim_size, second_dim_size],
                     order=order,
                 )
@@ -366,13 +379,13 @@ def initialise_3D_scikit_fem(self):
             u = self.u_sol[:, idx]
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, self.known_evals[t]
+                    t, u, self.inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u), [len_y, len_z]
+                    self.base_variable.evaluate(t, u, self.inputs), [len_y, len_z]
                 )
 
         # assign attributes for reference
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index e9260059eb..4ce57ce3bd 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -152,7 +152,7 @@ def build(self, check_model=True):
             self._model, inplace=False, check_model=check_model
         )
 
-    def solve(self, t_eval=None, solver=None, check_model=True):
+    def solve(self, t_eval=None, solver=None, inputs=None, check_model=True):
         """
         A method to solve the model. This method will automatically build
         and set the model parameters if not already done so.
@@ -166,6 +166,8 @@ def solve(self, t_eval=None, solver=None, check_model=True):
             non-dimensional time of 1.
         solver : :class:`pybamm.BaseSolver`
             The solver to use to solve the model.
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         check_model : bool, optional
             If True, model checks are performed after discretisation (see
             :meth:`pybamm.Discretisation.process_model`). Default is True.
@@ -185,9 +187,10 @@ def solve(self, t_eval=None, solver=None, check_model=True):
         if solver is None:
             solver = self.solver
 
-        self._solution = solver.solve(self.built_model, t_eval)
+        self.t_eval = t_eval
+        self._solution = solver.solve(self.built_model, t_eval, inputs=inputs)
 
-    def step(self, dt, solver=None, external_variables=None, save=True):
+    def step(self, dt, solver=None, external_variables=None, inputs=None, save=True):
         """
         A method to step the model forward one timestep. This method will
         automatically build and set the model parameters if not already done so.
@@ -203,6 +206,8 @@ def step(self, dt, solver=None, external_variables=None, save=True):
             values at the current time. The variables must correspond to
             the variables that would normally be found by solving the
             submodels that have been made external.
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         save : bool
             Turn on to store the solution of all previous timesteps
         """
@@ -212,7 +217,7 @@ def step(self, dt, solver=None, external_variables=None, save=True):
             solver = self.solver
 
         solution = solver.step(
-            self.built_model, dt, external_variables=external_variables
+            self.built_model, dt, external_variables=external_variables, inputs=inputs
         )
 
         if save is False or self._made_first_step is False:
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 5ef78b690e..7e93be8b58 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -49,7 +49,7 @@ def atol(self):
     def atol(self, value):
         self._atol = value
 
-    def solve(self, model, t_eval):
+    def solve(self, model, t_eval, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -61,6 +61,8 @@ def solve(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Raises
         ------
@@ -77,14 +79,17 @@ def solve(self, model, t_eval):
         # Set up
         timer = pybamm.Timer()
         start_time = timer.time()
+        inputs = inputs or {}
         if model.convert_to_format == "casadi" or isinstance(self, pybamm.CasadiSolver):
-            self.set_up_casadi(model)
+            self.set_up_casadi(model, inputs)
         else:
-            self.set_up(model)
+            self.set_up(model, inputs)
         set_up_time = timer.time() - start_time
 
         # Solve
-        solution, solve_time, termination = self.compute_solution(model, t_eval)
+        solution, solve_time, termination = self.compute_solution(
+            model, t_eval, inputs=inputs
+        )
 
         # Assign times
         solution.solve_time = solve_time
@@ -100,7 +105,7 @@ def solve(self, model, t_eval):
         )
         return solution
 
-    def step(self, model, dt, npts=2, log=True, external_variables=None):
+    def step(self, model, dt, npts=2, log=True, external_variables=None, inputs=None):
         """
         Step the solution of the model forward by a given time increment. The
         first time this method is called it executes the necessary setup by
@@ -119,6 +124,9 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
         external_variables : dict
             A dictionary of external variables and their corresponding
             values at the current time
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
 
         Raises
         ------
@@ -132,6 +140,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
 
         # Set timer
         timer = pybamm.Timer()
+        inputs = inputs or {}
 
         if not hasattr(self, "y0"):
             # create a y_pad vector of the correct size:
@@ -148,12 +157,12 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
             if model.convert_to_format == "casadi" or isinstance(
                 self, pybamm.CasadiSolver
             ):
-                self.set_up_casadi(model)
+                self.set_up_casadi(model, inputs)
             else:
                 pybamm.logger.debug(
                     "Start stepping {} with {}".format(model.name, self.name)
                 )
-                self.set_up(model)
+                self.set_up(model, inputs)
             self.t = 0.0
             set_up_time = timer.time()
 
@@ -162,7 +171,7 @@ def step(self, model, dt, npts=2, log=True, external_variables=None):
 
         # Step
         t_eval = np.linspace(self.t, self.t + dt, npts)
-        solution, solve_time, termination = self.compute_solution(model, t_eval)
+        solution, solve_time, termination = self.compute_solution(model, t_eval, inputs)
 
         # Set self.t and self.y0 to their values at the final step
         self.t = solution.t[-1]
@@ -221,7 +230,7 @@ def set_external_variables(self, model, external_variables):
                 )
             self.y_ext[y_slice] = var_vals
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times. Note: this
         does *not* execute the solver setup.
 
@@ -232,11 +241,13 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         raise NotImplementedError
 
-    def set_up(self, model):
+    def set_up(self, model, inputs=None):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -244,11 +255,13 @@ def set_up(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         raise NotImplementedError
 
-    def set_up_casadi(self, model):
+    def set_up_casadi(self, model, inputs=None):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -256,6 +269,8 @@ def set_up_casadi(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         raise NotImplementedError
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index b5c0b5c0c5..7c0329d622 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -68,7 +68,7 @@ def __init__(
         self.extra_options = extra_options
         self.name = "CasADi solver ({}) with '{}' mode".format(method, mode)
 
-    def solve(self, model, t_eval):
+    def solve(self, model, t_eval, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -80,6 +80,8 @@ def solve(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Raises
         ------
@@ -99,7 +101,8 @@ def solve(self, model, t_eval):
         elif self.mode == "safe":
             # Step-and-check
             timer = pybamm.Timer()
-            self.set_up_casadi(model)
+            inputs = inputs or {}
+            self.set_up_casadi(model, inputs)
             set_up_time = timer.time()
             init_event_signs = np.sign(
                 np.concatenate([event(0, self.y0) for event in self.event_funs])
@@ -175,7 +178,7 @@ def solve(self, model, t_eval):
             )
             return solution
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times. In this class, we
         overwrite the behaviour of :class:`pybamm.DaeSolver`, since CasADi requires
         slightly different syntax.
@@ -187,6 +190,8 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         timer = pybamm.Timer()
@@ -198,6 +203,7 @@ def compute_solution(self, model, t_eval):
             self.casadi_algebraic,
             self.y0,
             t_eval,
+            inputs,
             mass_matrix=model.mass_matrix.entries,
         )
         solve_time = timer.time() - solve_start_time
@@ -207,7 +213,9 @@ def compute_solution(self, model, t_eval):
 
         return solution, solve_time, termination
 
-    def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
+    def integrate_casadi(
+        self, rhs, algebraic, y0, t_eval, inputs=None, mass_matrix=None
+    ):
         """
         Solve a DAE model defined by residuals with initial conditions y0.
 
@@ -220,11 +228,14 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
             The initial conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         mass_matrix : array_like, optional
             The (sparse) mass matrix for the chosen spatial method. This is only passed
             to check that the mass matrix is diagonal with 1s for the odes and 0s for
             the algebraic equations, as CasADi does not allow to pass mass matrices.
         """
+        inputs = inputs or {}
         options = {
             "grid": t_eval,
             "reltol": self.rtol,
@@ -238,19 +249,15 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, mass_matrix=None):
 
         # set up and solve
         t = casadi.MX.sym("t")
-        y_diff = casadi.MX.sym("y_diff", rhs(0, y0).shape[0])
+        u = casadi.vertcat(*[x for x in inputs.values()])
+        y_diff = casadi.MX.sym("y_diff", rhs(0, y0, u).shape[0])
+        problem = {"t": t, "x": y_diff}
         if algebraic is None:
-            problem = {"t": t, "x": y_diff, "ode": rhs(t, y_diff)}
+            problem.update({"ode": rhs(t, y_diff, u)})
         else:
-            y_alg = casadi.MX.sym("y_alg", algebraic(0, y0).shape[0])
+            y_alg = casadi.MX.sym("y_alg", algebraic(0, y0, u).shape[0])
             y = casadi.vertcat(y_diff, y_alg)
-            problem = {
-                "t": t,
-                "x": y_diff,
-                "z": y_alg,
-                "ode": rhs(t, y),
-                "alg": algebraic(t, y),
-            }
+            problem.update({"z": y_alg, "ode": rhs(t, y, u), "alg": algebraic(t, y, u)})
         integrator = casadi.integrator("F", self.method, problem, options)
         try:
             # Try solving
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 33d312de49..1822074368 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -67,7 +67,7 @@ def max_steps(self):
     def max_steps(self, max_steps):
         self._max_steps = max_steps
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times.
 
         Parameters
@@ -77,18 +77,13 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
-
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
         """
         timer = pybamm.Timer()
 
-        # update y_pad and y_ext
-        self.rhs.set_pad_ext(self.y_pad, self.y_ext)
-        self.algebraic.set_pad_ext(self.y_pad, self.y_ext)
-        self.residuals.set_pad_ext(self.y_pad, self.y_ext)
-        for evnt in self.event_funs:
-            evnt.set_pad_ext(self.y_pad, self.y_ext)
-        if self.jacobian:
-            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+        # Set inputs and external
+        self.set_inputs_and_external(inputs)
 
         solve_start_time = timer.time()
         pybamm.logger.info("Calling DAE solver")
@@ -109,7 +104,7 @@ def compute_solution(self, model, t_eval):
 
         return solution, solve_time, termination
 
-    def set_up(self, model):
+    def set_up(self, model, inputs=None):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -117,7 +112,12 @@ def set_up(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
         """
+        inputs = inputs or {}
+
         # create simplified rhs, algebraic and event expressions
         concatenated_rhs = model.concatenated_rhs
         concatenated_algebraic = model.concatenated_algebraic
@@ -159,11 +159,9 @@ def set_up(self, model):
                 jac = pybamm.EvaluatorPython(jac)
 
             jacobian = Jacobian(jac.evaluate)
-
-            def jacobian_alg(t, y):
-                y = y[:, np.newaxis]
-                y = add_external(y, self.y_pad, self.y_ext)
-                return jac_algebraic.evaluate(t, y, known_evals={})[0]
+            jacobian_alg = JacobianAlgebraic(jac_algebraic.evaluate)
+            jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+            jacobian_alg.set_inputs(inputs)
 
         else:
             jacobian = None
@@ -183,6 +181,11 @@ def jacobian_alg(t, y):
         rhs = Rhs(concatenated_rhs.evaluate)
         algebraic = Algebraic(concatenated_algebraic.evaluate)
 
+        rhs.set_pad_ext(self.y_pad, self.y_ext)
+        rhs.set_inputs(inputs)
+        algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        algebraic.set_inputs(inputs)
+
         if len(model.algebraic) > 0:
             y0 = self.calculate_consistent_initial_conditions(
                 rhs,
@@ -209,7 +212,7 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in events.values()]
         self.jacobian = jacobian
 
-    def set_up_casadi(self, model):
+    def set_up_casadi(self, model, inputs=None):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -217,15 +220,22 @@ def set_up_casadi(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
         """
+        inputs = inputs or {}
+
         # Convert model attributes to casadi
         t_casadi = casadi.MX.sym("t")
         y0 = model.concatenated_initial_conditions
         y0 = add_external(y0, self.y_pad, self.y_ext)
 
-        y_diff = casadi.MX.sym("y_diff", len(model.concatenated_rhs.evaluate(0, y0)))
+        y_diff = casadi.MX.sym(
+            "y_diff", len(model.concatenated_rhs.evaluate(0, y0, inputs))
+        )
         y_alg = casadi.MX.sym(
-            "y_alg", len(model.concatenated_algebraic.evaluate(0, y0))
+            "y_alg", len(model.concatenated_algebraic.evaluate(0, y0, inputs))
         )
         y_casadi = casadi.vertcat(y_diff, y_alg)
         if self.y_pad is not None:
@@ -233,47 +243,55 @@ def set_up_casadi(self, model):
             y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
         else:
             y_casadi_w_ext = y_casadi
+        u_casadi = {name: casadi.MX.sym(name) for name in inputs.keys()}
 
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(
+            t_casadi, y_casadi_w_ext, u_casadi
+        )
         pybamm.logger.info("Converting algebraic to CasADi")
         concatenated_algebraic = model.concatenated_algebraic.to_casadi(
-            t_casadi, y_casadi_w_ext
+            t_casadi, y_casadi_w_ext, u_casadi
         )
         all_states = casadi.vertcat(concatenated_rhs, concatenated_algebraic)
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi_w_ext)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext, u_casadi)
             for name, event in model.events.items()
         }
 
         # Create functions to evaluate rhs and algebraic
+        u_casadi_stacked = casadi.vertcat(*[u for u in u_casadi.values()])
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [concatenated_rhs]
         )
         concatenated_algebraic_fn = casadi.Function(
-            "algebraic", [t_casadi, y_casadi_w_ext], [concatenated_algebraic]
+            "algebraic",
+            [t_casadi, y_casadi_w_ext, u_casadi_stacked],
+            [concatenated_algebraic],
+        )
+        all_states_fn = casadi.Function(
+            "all", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [all_states]
         )
-        all_states_fn = casadi.Function("all", [t_casadi, y_casadi_w_ext], [all_states])
 
         if model.use_jacobian:
 
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(all_states, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [casadi_jac]
             )
             casadi_jac_alg = casadi.jacobian(concatenated_algebraic, y_casadi)
             casadi_jac_alg_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac_alg]
+                "jacobian",
+                [t_casadi, y_casadi_w_ext, u_casadi_stacked],
+                [casadi_jac_alg],
             )
 
             jacobian = JacobianCasadi(casadi_jac_fn)
-
-            def jacobian_alg(t, y):
-                y = y[:, np.newaxis]
-                y = add_external(y, self.y_pad, self.y_ext)
-                return casadi_jac_alg_fn(t, y)
+            jacobian_alg = JacobianAlgebraicCasadi(casadi_jac_alg_fn)
+            jacobian_alg.set_pad_ext(self.y_pad, self.y_ext)
+            jacobian_alg.set_inputs(inputs)
 
         else:
             jacobian = None
@@ -283,7 +301,9 @@ def jacobian_alg(t, y):
         algebraic = AlgebraicCasadi(concatenated_algebraic_fn)
 
         rhs.set_pad_ext(self.y_pad, self.y_ext)
+        rhs.set_inputs(inputs)
         algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        algebraic.set_inputs(inputs)
 
         if len(model.algebraic) > 0:
 
@@ -300,9 +320,9 @@ def jacobian_alg(t, y):
         # Create event-dependent function to evaluate events
         def get_event_class(event):
             casadi_event_fn = casadi.Function(
-                "event", [t_casadi, y_casadi_w_ext], [event]
+                "event", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [event]
             )
-            return EvalEvent(casadi_event_fn)
+            return EvalEventCasadi(casadi_event_fn)
 
         # Add the solver attributes
         # Note: these are the (possibly) converted to python version rhs, algebraic
@@ -319,6 +339,30 @@ def get_event_class(event):
         self.casadi_rhs = concatenated_rhs_fn
         self.casadi_algebraic = concatenated_algebraic_fn
 
+    def set_inputs_and_external(self, inputs):
+        """
+        Set values that are controlled externally, such as external variables and input
+        parameters
+
+        Parameters
+        ----------
+        inputs : dict
+            Any input parameters to pass to the model when solving
+
+        """
+        self.rhs.set_pad_ext(self.y_pad, self.y_ext)
+        self.rhs.set_inputs(inputs)
+        self.algebraic.set_pad_ext(self.y_pad, self.y_ext)
+        self.algebraic.set_inputs(inputs)
+        self.residuals.set_pad_ext(self.y_pad, self.y_ext)
+        self.residuals.set_inputs(inputs)
+        for evnt in self.event_funs:
+            evnt.set_pad_ext(self.y_pad, self.y_ext)
+            evnt.set_inputs(inputs)
+        if self.jacobian:
+            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+            self.jacobian.set_inputs(inputs)
+
     def calculate_consistent_initial_conditions(
         self, rhs, algebraic, y0_guess, jac=None
     ):
@@ -439,22 +483,32 @@ def integrate(
         raise NotImplementedError
 
 
-class Rhs:
-    "Returns information about rhs at time t and state y"
-
-    def __init__(self, concatenated_rhs_fn):
-        self.concatenated_rhs_fn = concatenated_rhs_fn
-        self.y_pad = None
-        self.y_ext = None
+class SolverCallable:
+    "A class that will be called by the solver when integrating"
+    y_pad = None
+    y_ext = None
+    inputs = {}
+    inputs_casadi = casadi.DM()
 
     def set_pad_ext(self, y_pad, y_ext):
         self.y_pad = y_pad
         self.y_ext = y_ext
 
+    def set_inputs(self, inputs):
+        self.inputs = inputs
+        self.inputs_casadi = casadi.vertcat(*[x for x in inputs.values()])
+
+
+class Rhs(SolverCallable):
+    "Returns information about rhs at time t and state y"
+
+    def __init__(self, concatenated_rhs_fn):
+        self.concatenated_rhs_fn = concatenated_rhs_fn
+
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_rhs_fn(t, y, known_evals={})[0][:, 0]
+        return self.concatenated_rhs_fn(t, y, self.inputs, known_evals={})[0][:, 0]
 
 
 class RhsCasadi(Rhs):
@@ -463,25 +517,21 @@ class RhsCasadi(Rhs):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_rhs_fn(t, y).full()[:, 0]
+        return self.concatenated_rhs_fn(t, y, self.inputs_casadi).full()[:, 0]
 
 
-class Algebraic:
+class Algebraic(SolverCallable):
     "Returns information about algebraic equations at time t and state y"
 
     def __init__(self, concatenated_algebraic_fn):
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_algebraic_fn(t, y, known_evals={})[0][:, 0]
+        return self.concatenated_algebraic_fn(t, y, self.inputs, known_evals={})[0][
+            :, 0
+        ]
 
 
 class AlgebraicCasadi(Algebraic):
@@ -490,10 +540,10 @@ class AlgebraicCasadi(Algebraic):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.concatenated_algebraic_fn(t, y).full()[:, 0]
+        return self.concatenated_algebraic_fn(t, y, self.inputs_casadi).full()[:, 0]
 
 
-class Residuals:
+class Residuals(SolverCallable):
     "Returns information about residuals at time t and state y"
 
     def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn):
@@ -501,12 +551,6 @@ def __init__(self, model, concatenated_rhs_fn, concatenated_algebraic_fn):
         self.concatenated_rhs_fn = concatenated_rhs_fn
         self.concatenated_algebraic_fn = concatenated_algebraic_fn
         self.mass_matrix = model.mass_matrix.entries
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y, ydot):
         pybamm.logger.debug(
@@ -514,9 +558,13 @@ def __call__(self, t, y, ydot):
         )
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        rhs_eval, known_evals = self.concatenated_rhs_fn(t, y, known_evals={})
+        rhs_eval, known_evals = self.concatenated_rhs_fn(
+            t, y, self.inputs, known_evals={}
+        )
         # reuse known_evals
-        alg_eval = self.concatenated_algebraic_fn(t, y, known_evals=known_evals)[0]
+        alg_eval = self.concatenated_algebraic_fn(
+            t, y, self.inputs, known_evals=known_evals
+        )[0]
         # turn into 1D arrays
         rhs_eval = rhs_eval[:, 0]
         alg_eval = alg_eval[:, 0]
@@ -530,8 +578,6 @@ def __init__(self, model, all_states_fn):
         self.model = model
         self.all_states_fn = all_states_fn
         self.mass_matrix = model.mass_matrix.entries
-        self.y_pad = None
-        self.y_ext = None
 
     def __call__(self, t, y, ydot):
         pybamm.logger.debug(
@@ -539,44 +585,44 @@ def __call__(self, t, y, ydot):
         )
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        states_eval = self.all_states_fn(t, y).full()[:, 0]
+        states_eval = self.all_states_fn(t, y, self.inputs_casadi).full()[:, 0]
         return states_eval - self.mass_matrix @ ydot
 
 
-class EvalEvent:
+class EvalEvent(SolverCallable):
     "Returns information about events at time t and state y"
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
-        self.y_pad = None
-        self.y_ext = None
 
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.event_fn(t, y, self.inputs)
+
+
+class EvalEventCasadi(EvalEvent):
+    "Returns information about events at time t and state y"
+
+    def __init__(self, event_fn):
+        self.event_fn = event_fn
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.event_fn(t, y)
+        return self.event_fn(t, y, self.inputs_casadi)
 
 
-class Jacobian:
+class Jacobian(SolverCallable):
     "Returns information about the jacobian at time t and state y"
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y, known_evals={})[0]
+        return self.jac_fn(t, y, self.inputs, known_evals={})[0]
 
 
 class JacobianCasadi(Jacobian):
@@ -585,5 +631,28 @@ class JacobianCasadi(Jacobian):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y)
+        return self.jac_fn(t, y, self.inputs_casadi)
+
+
+class JacobianAlgebraic(SolverCallable):
+    "Returns information about the algebraic part of the jacobian at time t and state y"
 
+    def __init__(self, jac_alg_fn):
+        self.jac_fn = jac_alg_fn
+
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.jac_fn(t, y, self.inputs, known_evals={})[0]
+
+
+class JacobianAlgebraicCasadi(JacobianAlgebraic):
+    """
+    Returns information about the algebraic part of the jacobian at time t and
+    state y, with CasADi
+    """
+
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.jac_fn(t, y, self.inputs_casadi)
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 934e521ef5..8faf202cdc 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -23,7 +23,7 @@ def __init__(self, method=None, rtol=1e-6, atol=1e-6):
         super().__init__(method, rtol, atol)
         self.name = "Base ODE solver"
 
-    def compute_solution(self, model, t_eval):
+    def compute_solution(self, model, t_eval, inputs=None):
         """Calculate the solution of the model at specified times.
 
         Parameters
@@ -33,16 +33,16 @@ def compute_solution(self, model, t_eval):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         """
         timer = pybamm.Timer()
 
-        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
-        for evnt in self.event_funs:
-            evnt.set_pad_ext(self.y_pad, self.y_ext)
-        if self.jacobian:
-            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+        # Set inputs and external
+        self.set_inputs_and_external(inputs)
 
+        # Solve
         solve_start_time = timer.time()
         pybamm.logger.info("Calling ODE solver")
         solution = self.integrate(
@@ -61,7 +61,7 @@ def compute_solution(self, model, t_eval):
 
         return solution, solve_time, termination
 
-    def set_up(self, model):
+    def set_up(self, model, inputs=None):
         """Unpack model, perform checks, simplify and calculate jacobian.
 
         Parameters
@@ -69,6 +69,9 @@ def set_up(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
+
 
         Raises
         ------
@@ -83,6 +86,8 @@ def set_up(self, model):
                 """Cannot use ODE solver to solve model with DAEs"""
             )
 
+        inputs = inputs or {}
+
         # create simplified rhs and event expressions
         concatenated_rhs = model.concatenated_rhs
         events = model.events
@@ -146,7 +151,7 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in events.values()]
         self.jacobian = jacobian
 
-    def set_up_casadi(self, model):
+    def set_up_casadi(self, model, inputs=None):
         """Convert model to casadi format and use their inbuilt functionalities.
 
         Parameters
@@ -154,6 +159,8 @@ def set_up_casadi(self, model):
         model : :class:`pybamm.BaseModel`
             The model whose solution to calculate. Must have attributes rhs and
             initial_conditions
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Raises
         ------
@@ -172,6 +179,8 @@ def set_up_casadi(self, model):
 
         t_casadi = casadi.MX.sym("t")
         y_casadi = casadi.MX.sym("y", len(y0))
+        inputs = inputs or {}
+        u_casadi = {name: casadi.MX.sym(name) for name in inputs.keys()}
 
         if self.y_pad is not None:
             y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
@@ -180,31 +189,34 @@ def set_up_casadi(self, model):
             y_casadi_w_ext = y_casadi
 
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(
+            t_casadi, y_casadi_w_ext, u_casadi
+        )
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi_w_ext)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext, u_casadi)
             for name, event in model.events.items()
         }
 
         # Create function to evaluate rhs
+        u_casadi_stacked = casadi.vertcat(*[u for u in u_casadi.values()])
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [concatenated_rhs]
         )
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
             casadi_event_fn = casadi.Function(
-                "event", [t_casadi, y_casadi_w_ext], [event]
+                "event", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [event]
             )
-            return EvalEvent(casadi_event_fn)
+            return EvalEventCasadi(casadi_event_fn)
 
         # Create function to evaluate jacobian
         if model.use_jacobian:
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(concatenated_rhs, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext, u_casadi_stacked], [casadi_jac]
             )
 
             jacobian = JacobianCasadi(casadi_jac_fn)
@@ -219,6 +231,26 @@ def get_event_class(event):
         self.event_funs = [get_event_class(event) for event in casadi_events.values()]
         self.jacobian = jacobian
 
+    def set_inputs_and_external(self, inputs):
+        """
+        Set values that are controlled externally, such as external variables and input
+        parameters
+
+        Parameters
+        ----------
+        inputs : dict
+            Any input parameters to pass to the model when solving
+
+        """
+        self.dydt.set_pad_ext(self.y_pad, self.y_ext)
+        self.dydt.set_inputs(inputs)
+        for evnt in self.event_funs:
+            evnt.set_pad_ext(self.y_pad, self.y_ext)
+            evnt.set_inputs(inputs)
+        if self.jacobian:
+            self.jacobian.set_pad_ext(self.y_pad, self.y_ext)
+            self.jacobian.set_inputs(inputs)
+
     def integrate(
         self, derivs, y0, t_eval, events=None, mass_matrix=None, jacobian=None
     ):
@@ -244,25 +276,35 @@ def integrate(
         raise NotImplementedError
 
 
+class SolverCallable:
+    "A class that will be called by the solver when integrating"
+    y_pad = None
+    y_ext = None
+    inputs = {}
+    inputs_casadi = casadi.DM()
+
+    def set_pad_ext(self, y_pad, y_ext):
+        self.y_pad = y_pad
+        self.y_ext = y_ext
+
+    def set_inputs(self, inputs):
+        self.inputs = inputs
+        self.inputs_casadi = casadi.vertcat(*[x for x in inputs.values()])
+
+
 # Set up caller classes outside of the solver object to allow pickling
-class Dydt:
+class Dydt(SolverCallable):
     "Returns information about time derivatives at time t and state y"
 
     def __init__(self, model, concatenated_rhs_fn):
         self.model = model
         self.concatenated_rhs_fn = concatenated_rhs_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         pybamm.logger.debug("Evaluating RHS for {} at t={}".format(self.model.name, t))
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        dy = self.concatenated_rhs_fn(t, y, known_evals={})[0]
+        dy = self.concatenated_rhs_fn(t, y, self.inputs, known_evals={})[0]
         return dy[:, 0]
 
 
@@ -273,44 +315,44 @@ def __call__(self, t, y):
         pybamm.logger.debug("Evaluating RHS for {} at t={}".format(self.model.name, t))
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        dy = self.concatenated_rhs_fn(t, y).full()
+        dy = self.concatenated_rhs_fn(t, y, self.inputs_casadi).full()
         return dy[:, 0]
 
 
-class EvalEvent:
+class EvalEvent(SolverCallable):
     "Returns information about events at time t and state y"
 
     def __init__(self, event_fn):
         self.event_fn = event_fn
-        self.y_pad = None
-        self.y_ext = None
 
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
+    def __call__(self, t, y):
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
+        return self.event_fn(t, y, self.inputs)
+
+
+class EvalEventCasadi(EvalEvent):
+    "Returns information about events at time t and state y"
+
+    def __init__(self, event_fn):
+        self.event_fn = event_fn
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.event_fn(t, y)
+        return self.event_fn(t, y, self.inputs_casadi)
 
 
-class Jacobian:
+class Jacobian(SolverCallable):
     "Returns information about the jacobian at time t and state y"
 
     def __init__(self, jac_fn):
         self.jac_fn = jac_fn
-        self.y_pad = None
-        self.y_ext = None
-
-    def set_pad_ext(self, y_pad, y_ext):
-        self.y_pad = y_pad
-        self.y_ext = y_ext
 
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y, known_evals={})[0]
+        return self.jac_fn(t, y, self.inputs, known_evals={})[0]
 
 
 class JacobianCasadi(Jacobian):
@@ -319,4 +361,4 @@ class JacobianCasadi(Jacobian):
     def __call__(self, t, y):
         y = y[:, np.newaxis]
         y = add_external(y, self.y_pad, self.y_ext)
-        return self.jac_fn(t, y)
+        return self.jac_fn(t, y, self.inputs_casadi)
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index fd2af8d763..1541a3b801 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -47,6 +47,7 @@ def integrate(
         jacobian : method, optional
             A function that takes in t and y and returns the Jacobian. If
             None, the solver will approximate the Jacobian.
+
         Returns
         -------
         object
diff --git a/tests/unit/test_expression_tree/test_input_parameter.py b/tests/unit/test_expression_tree/test_input_parameter.py
new file mode 100644
index 0000000000..f5c24e3909
--- /dev/null
+++ b/tests/unit/test_expression_tree/test_input_parameter.py
@@ -0,0 +1,38 @@
+#
+# Tests for the InputParameter class
+#
+import numbers
+import pybamm
+import unittest
+
+
+class TestInputParameter(unittest.TestCase):
+    def test_input_parameter_init(self):
+        a = pybamm.InputParameter("a")
+        self.assertEqual(a.name, "a")
+        self.assertEqual(a.evaluate(u={"a": 1}), 1)
+        self.assertEqual(a.evaluate(u={"a": 5}), 5)
+
+    def test_evaluate_for_shape(self):
+        a = pybamm.InputParameter("a")
+        self.assertIsInstance(a.evaluate_for_shape(), numbers.Number)
+
+    def test_errors(self):
+        a = pybamm.InputParameter("a")
+        with self.assertRaises(TypeError):
+            a.evaluate(u="not a dictionary")
+        with self.assertRaises(KeyError):
+            a.evaluate(u={"bad param": 5})
+        # if u is not provided it gets turned into a dictionary and then raises KeyError
+        with self.assertRaises(KeyError):
+            a.evaluate()
+
+
+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_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index ab7d50d9b9..d0ee4fcf1e 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -157,6 +157,40 @@ def myfunction(x, y):
         f = pybamm.Function(myfunction, a, b).diff(b)
         self.assert_casadi_equal(f.to_casadi(), casadi.MX(3), evalf=True)
 
+    def test_convert_input_parameter(self):
+        # Arrays
+        a = np.array([1, 2, 3, 4, 5])
+        pybamm_a = pybamm.Array(a)
+        self.assert_casadi_equal(pybamm_a.to_casadi(), casadi.MX(a))
+
+        casadi_t = casadi.MX.sym("t")
+        casadi_y = casadi.MX.sym("y", 10)
+        casadi_us = {
+            "Input 1": casadi.MX.sym("Input 1"),
+            "Input 2": casadi.MX.sym("Input 2"),
+        }
+
+        pybamm_y = pybamm.StateVector(slice(0, 10))
+        pybamm_u1 = pybamm.InputParameter("Input 1")
+        pybamm_u2 = pybamm.InputParameter("Input 2")
+
+        # Input only
+        self.assert_casadi_equal(
+            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_us), casadi_us["Input 1"]
+        )
+
+        # More complex
+        expr = pybamm_u1 + pybamm_y
+        self.assert_casadi_equal(
+            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            casadi_us["Input 1"] + casadi_y,
+        )
+        expr = pybamm_u2 * pybamm_y
+        self.assert_casadi_equal(
+            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            casadi_us["Input 2"] * casadi_y,
+        )
+
     def test_errors(self):
         y = pybamm.StateVector(slice(0, 10))
         with self.assertRaisesRegex(
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index 3161967d20..c969e89928 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -37,6 +37,7 @@ def test_symbol_new_copy(self):
             pybamm.NumpyConcatenation(a, b, v_s),
             pybamm.DomainConcatenation([v_n, v_s], mesh),
             pybamm.Parameter("param"),
+            pybamm.InputParameter("param"),
             pybamm.StateVector(slice(0, 56)),
             pybamm.Matrix(np.ones((50, 40))),
             pybamm.SpatialVariable("x", ["negative electrode"]),
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index 05cd23b7da..c5182c3c8c 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -228,6 +228,23 @@ def test_process_symbol(self):
         with self.assertRaises(NotImplementedError):
             parameter_values.process_symbol(sym)
 
+    def test_process_input_parameter(self):
+        parameter_values = pybamm.ParameterValues({"a": "[input]", "b": 3})
+        # process input parameter
+        a = pybamm.Parameter("a")
+        processed_a = parameter_values.process_symbol(a)
+        self.assertIsInstance(processed_a, pybamm.InputParameter)
+        self.assertEqual(processed_a.evaluate(u={"a": 5}), 5)
+
+        # process binary operation
+        b = pybamm.Parameter("b")
+        add = a + b
+        processed_add = parameter_values.process_symbol(add)
+        self.assertIsInstance(processed_add, pybamm.Addition)
+        self.assertIsInstance(processed_add.children[0], pybamm.InputParameter)
+        self.assertIsInstance(processed_add.children[1], pybamm.Scalar)
+        self.assertEqual(processed_add.evaluate(u={"a": 4}), 7)
+
     def test_process_function_parameter(self):
         parameter_values = pybamm.ParameterValues(
             {
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 13d66ff0d2..74cee9d218 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -16,8 +16,9 @@ def test_integrate(self):
 
         t = casadi.MX.sym("t")
         y = casadi.MX.sym("y")
+        u = casadi.MX.sym("u")
         constant_growth = casadi.MX(0.5)
-        rhs = casadi.Function("rhs", [t, y], [constant_growth])
+        rhs = casadi.Function("rhs", [t, y, u], [constant_growth])
 
         y0 = np.array([0])
         t_eval = np.linspace(0, 1, 100)
@@ -29,7 +30,7 @@ def test_integrate(self):
         solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="cvodes")
 
         exponential_decay = -0.1 * y
-        rhs = casadi.Function("rhs", [t, y], [exponential_decay])
+        rhs = casadi.Function("rhs", [t, y, u], [exponential_decay])
 
         y0 = np.array([1])
         t_eval = np.linspace(0, 1, 100)
@@ -37,18 +38,31 @@ def test_integrate(self):
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
         self.assertEqual(solution.termination, "final time")
 
+        # Exponential decay with input
+        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="cvodes")
+
+        exponential_decay = -u * y
+        rhs = casadi.Function("rhs", [t, y, u], [exponential_decay])
+
+        y0 = np.array([1])
+        t_eval = np.linspace(0, 1, 100)
+        solution = solver.integrate_casadi(rhs, None, y0, t_eval, inputs={"u": 0.1})
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+        self.assertEqual(solution.termination, "final time")
+
     def test_integrate_failure(self):
         # Turn off warnings to ignore sqrt error
         warnings.simplefilter("ignore")
 
         t = casadi.MX.sym("t")
         y = casadi.MX.sym("y")
+        u = casadi.MX.sym("u")
         sqrt_decay = -np.sqrt(y)
 
         y0 = np.array([1])
         t_eval = np.linspace(0, 3, 100)
         solver = pybamm.CasadiSolver()
-        rhs = casadi.Function("rhs", [t, y], [sqrt_decay])
+        rhs = casadi.Function("rhs", [t, y, u], [sqrt_decay])
         # Expect solver to fail when y goes negative
         with self.assertRaises(pybamm.SolverError):
             solver.integrate_casadi(rhs, None, y0, t_eval)
@@ -184,6 +198,30 @@ def test_model_step(self):
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], step_sol.y[0])
 
+    def test_model_solver_with_inputs(self):
+        # Create model
+        model = pybamm.BaseModel()
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = {"var=0.5": pybamm.min(var - 0.5)}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
+        self.assertLess(len(solution.t), len(t_eval))
+        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_dae_solver.py b/tests/unit/test_solvers/test_dae_solver.py
index 6e35964415..4eb334e722 100644
--- a/tests/unit/test_solvers/test_dae_solver.py
+++ b/tests/unit/test_solvers/test_dae_solver.py
@@ -77,6 +77,11 @@ def algebraic(t, y):
         ):
             solver.calculate_consistent_initial_conditions(rhs, algebraic, y0)
 
+    def test_errors(self):
+        solver = pybamm.DaeSolver()
+        with self.assertRaises(NotImplementedError):
+            solver.integrate(None, None, None)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 7ac16b5af3..c2259b6a79 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -576,6 +576,7 @@ def jacobian(t, y):
                 ]
             )
 
+        model.jacobian = jacobian
         # Solve
         solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
         t_eval = np.linspace(0, 1, 100)
@@ -717,6 +718,33 @@ def test_model_solver_dae_events_casadi(self):
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
             np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
+    def test_model_solver_dae_inputs_events(self):
+        # Create model
+        for form in ["python", "casadi"]:
+            model = pybamm.BaseModel()
+            model.convert_to_format = form
+            whole_cell = ["negative electrode", "separator", "positive electrode"]
+            var1 = pybamm.Variable("var1", domain=whole_cell)
+            var2 = pybamm.Variable("var2", domain=whole_cell)
+            model.rhs = {var1: pybamm.InputParameter("rate 1") * var1}
+            model.algebraic = {var2: pybamm.InputParameter("rate 2") * var1 - var2}
+            model.initial_conditions = {var1: 1, var2: 2}
+            model.events = {
+                "var1 = 1.5": pybamm.min(var1 - 1.5),
+                "var2 = 2.5": pybamm.min(var2 - 2.5),
+            }
+            disc = get_discretisation_for_testing()
+            disc.process_model(model)
+
+            # Solve
+            solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+            t_eval = np.linspace(0, 5, 100)
+            solution = solver.solve(model, t_eval, inputs={"rate 1": 0.1, "rate 2": 2})
+            np.testing.assert_array_less(solution.y[0], 1.5)
+            np.testing.assert_array_less(solution.y[-1], 2.5)
+            np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
+            np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
+
     def test_solve_ode_model_with_dae_solver_casadi(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "casadi"
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 68ecd1de05..ff10eadc8f 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -272,6 +272,31 @@ def test_model_step_python(self):
         solution = solver.solve(model, t_eval)
         np.testing.assert_allclose(solution.y[0], step_sol.y[0])
 
+    def test_model_solver_with_inputs(self):
+        # Create model
+        model = pybamm.BaseModel()
+        model.convert_to_format = "python"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = {"var=0.5": pybamm.min(var - 0.5)}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
+        self.assertLess(len(solution.t), len(t_eval))
+        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+
     def test_model_solver_with_event_with_casadi(self):
         # Create model
         model = pybamm.BaseModel()
@@ -299,6 +324,31 @@ def test_model_solver_with_event_with_casadi(self):
             np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
             np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
 
+    def test_model_solver_with_inputs_with_casadi(self):
+        # Create model
+        model = pybamm.BaseModel()
+        model.convert_to_format = "casadi"
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var = pybamm.Variable("var", domain=domain)
+        model.rhs = {var: -pybamm.InputParameter("rate") * var}
+        model.initial_conditions = {var: 1}
+        model.events = {"var=0.5": pybamm.min(var - 0.5)}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, inputs={"rate": 0.1})
+        self.assertLess(len(solution.t), len(t_eval))
+        np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
+        np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")