Repository: pybamm-team/PyBaMM

Issue:
Input external temperature matrix 
**Summary**
Allow for an external temperature matrix T(x,y,z) (and a potetial boundary conditions phi(y,z)) and output an internal heat generation matrix. 



Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index aac622f6d..bcaeaa1f0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,11 +2,14 @@
 
 ## Features
 
+-   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))
 -   Simplified the interface for setting current functions ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))
 -   Added Heaviside operator ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))
+-   New extrapolation methods ([#707](https://github.com/pybamm-team/PyBaMM/pull/707))
 -   Added some "Getting Started" documentation ([#703](https://github.com/pybamm-team/PyBaMM/pull/703))
+-   Allow abs tolerance to be set by variable for IDA KLU solver ([#700](https://github.com/pybamm-team/PyBaMM/pull/700))
 -   Added Simulation class ([#693](https://github.com/pybamm-team/PyBaMM/pull/693)) with load/save functionality ([#732](https://github.com/pybamm-team/PyBaMM/pull/732))
 -   Added interface to CasADi solver ([#687](https://github.com/pybamm-team/PyBaMM/pull/687), [#691](https://github.com/pybamm-team/PyBaMM/pull/691), [#714](https://github.com/pybamm-team/PyBaMM/pull/714)). This makes the SUNDIALS DAE solvers (Scikits and KLU) truly optional (though IDA KLU is recommended for solving the DFN).
 -   Added option to use CasADi's Algorithmic Differentiation framework to calculate Jacobians ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 690be6e08..23f55aa76 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -120,6 +120,7 @@ from .expression_tree.exceptions import (
     ModelWarning,
     UndefinedOperationError,
     GeometryError,
+    InputError,
 )
 
 # Operations
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 708fbeb50..13772f0c6 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -52,6 +52,7 @@ class Discretisation(object):
         self.bcs = {}
         self.y_slices = {}
         self._discretised_symbols = {}
+        self.external_variables = []
 
     @property
     def mesh(self):
@@ -119,12 +120,20 @@ class Discretisation(object):
 
         # Prepare discretisation
         # set variables (we require the full variable not just id)
-        variables = list(model.rhs.keys()) + list(model.algebraic.keys())
+        variables = (
+            list(model.rhs.keys())
+            + list(model.algebraic.keys())
+            + model.external_variables
+        )
 
         # Set the y split for variables
         pybamm.logger.info("Set variable slices for {}".format(model.name))
         self.set_variable_slices(variables)
 
+        # now add extrapolated external variables to the boundary conditions
+        # if required by the spatial method
+        self._preprocess_external_variables(model)
+
         # set boundary conditions (only need key ids for boundary_conditions)
         pybamm.logger.info("Discretise boundary conditions for {}".format(model.name))
         self.bcs = self.process_boundary_conditions(model)
@@ -142,7 +151,28 @@ class Discretisation(object):
 
         model_disc.bcs = self.bcs
 
-        # Process initial condtions
+        self.external_variables = model.external_variables
+        # find where external variables begin in state vector
+        # we always append external variables to the end, so
+        # it is sufficient to only know the starting location
+        start_vals = []
+        for var in self.external_variables:
+            if isinstance(var, pybamm.Concatenation):
+                for child in var.children:
+                    start_vals += [self.y_slices[child.id][0].start]
+            elif isinstance(var, pybamm.Variable):
+                start_vals += [self.y_slices[var.id][0].start]
+
+        # attach properties of the state vector so that it
+        # can be divided correctly during the solving stage
+        model_disc.external_variables = model.external_variables
+        model_disc.y_length = self.y_length
+        model_disc.y_slices = self.y_slices
+        if start_vals:
+            model_disc.external_start = min(start_vals)
+        else:
+            model_disc.external_start = self.y_length
+
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
         ics, concat_ics = self.process_initial_conditions(model)
         model_disc.initial_conditions = ics
@@ -223,10 +253,33 @@ class Discretisation(object):
                 start = end
 
         self.y_slices = y_slices
+        self.y_length = end
 
         # reset discretised_symbols
         self._discretised_symbols = {}
 
+    def _preprocess_external_variables(self, model):
+        """
+        A method to preprocess external variables so that they are
+        compatible with the spatial method. For example, in finite
+        volume, the user will supply a vector of values valid on the
+        cell centres. However, for model processing, we also require
+        the boundary edge fluxes. Therefore, we extrapolate and add
+        the boundary fluxes to the boundary conditions, which are
+        employed in generating the grad and div matrices.
+        The processing is delegated to spatial methods as
+        the preprocessing required for finite volume and finite
+        element will be different.
+        """
+
+        for var in model.external_variables:
+            if var.domain != []:
+                new_bcs = self.spatial_methods[
+                    var.domain[0]
+                ].preprocess_external_variables(var)
+
+                model.boundary_conditions.update(new_bcs)
+
     def set_internal_boundary_conditions(self, model):
         """
         A method to set the internal boundary conditions for the submodel.
@@ -814,7 +867,14 @@ class Discretisation(object):
         if check_complete:
             # Check keys from the given var_eqn_dict against self.y_slices
             ids = {v.id for v in unpacked_variables}
-            if ids != self.y_slices.keys():
+            external_id = {v.id for v in self.external_variables}
+            for var in self.external_variables:
+                child_ids = {child.id for child in var.children}
+                external_id = external_id.union(child_ids)
+            y_slices_with_external_removed = set(self.y_slices.keys()).difference(
+                external_id
+            )
+            if ids != y_slices_with_external_removed:
                 given_variable_names = [v.name for v in var_eqn_dict.keys()]
                 raise pybamm.ModelError(
                     "Initial conditions are insufficient. Only "
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index f00720134..cbeca72b5 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -59,3 +59,11 @@ class UndefinedOperationError(Exception):
     """
 
     pass
+
+
+class InputError(Exception):
+    """
+    An external variable has been input incorrectly into PyBaMM
+    """
+
+    pass
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index ecf51e994..78625704d 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -109,6 +109,7 @@ class BaseModel(object):
         self._mass_matrix = None
         self._jacobian = None
         self._jacobian_algebraic = None
+        self.external_variables = []
 
         # Default behaviour is to use the jacobian and simplify
         self.use_jacobian = True
@@ -400,7 +401,17 @@ class BaseModel(object):
         # If any variables in the equations don't appear in the keys then the model is
         # underdetermined
         vars_in_keys = vars_in_rhs_keys.union(vars_in_algebraic_keys)
-        extra_variables = vars_in_eqns.difference(vars_in_keys)
+        extra_variables_in_equations = vars_in_eqns.difference(vars_in_keys)
+
+        # get ids of external variables
+        external_ids = {var.id for var in self.external_variables}
+        for var in self.external_variables:
+            if isinstance(var, pybamm.Concatenation):
+                child_ids = {child.id for child in var.children}
+                external_ids = external_ids.union(child_ids)
+
+        extra_variables = extra_variables_in_equations.difference(external_ids)
+
         if extra_variables:
             raise pybamm.ModelError("model is underdetermined (too many variables)")
 
@@ -491,19 +502,27 @@ class BaseModel(object):
                 {x.id: x for x in eqn.pre_order() if isinstance(x, pybamm.Variable)}
             )
         var_ids_in_keys = set()
-        for var in {**self.rhs, **self.algebraic}.keys():
+
+        model_and_external_variables = (
+            list(self.rhs.keys())
+            + list(self.algebraic.keys())
+            + self.external_variables
+        )
+
+        for var in model_and_external_variables:
             if isinstance(var, pybamm.Variable):
                 var_ids_in_keys.add(var.id)
             # Key can be a concatenation
             elif isinstance(var, pybamm.Concatenation):
                 var_ids_in_keys.update([child.id for child in var.children])
+
         for var_id, var in all_vars.items():
             if var_id not in var_ids_in_keys:
                 raise pybamm.ModelError(
                     """
                     No key set for variable '{}'. Make sure it is included in either
-                    model.rhs or model.algebraic in an unmodified form (e.g. not
-                    Broadcasted)
+                    model.rhs, model.algebraic, or model.external_variables in an
+                    unmodified form (e.g. not Broadcasted)
                     """.format(
                         var
                     )
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 90a381447..b16e88a63 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -57,6 +57,11 @@ class BaseBatteryModel(pybamm.BaseModel):
                 only takes effect if "dimensionality" is 0. If "dimensionality"
                 is 1 or 2 current collector effects are always included. Must be 'False'
                 for lead-acid models.
+            * "external submodels" : list
+                A list of the submodels that you would like to supply an external
+                variable for instead of solving in PyBaMM. The entries of the lists
+                are strings that correspond to the submodel names in the keys
+                of `self.submodels`.
 
 
     **Extends:** :class:`pybamm.BaseModel`
@@ -68,6 +73,7 @@ class BaseBatteryModel(pybamm.BaseModel):
         self.set_standard_output_variables()
         self.submodels = {}
         self._built = False
+        self._built_fundamental_and_external = False
 
     @property
     def default_parameter_values(self):
@@ -150,6 +156,7 @@ class BaseBatteryModel(pybamm.BaseModel):
             "particle": "Fickian diffusion",
             "thermal": "isothermal",
             "thermal current collector": False,
+            "external submodels": []
         }
         options = default_options
         # any extra options overwrite the default options
@@ -389,17 +396,7 @@ class BaseBatteryModel(pybamm.BaseModel):
                 {"y": var.y, "y [m]": var.y * L_y, "z": var.z, "z [m]": var.z * L_z}
             )
 
-    def build_model(self):
-
-        # Check if already built
-        if self._built:
-            raise pybamm.ModelError(
-                """Model already built. If you are adding a new submodel, try using
-                `model.update` instead."""
-            )
-
-        pybamm.logger.info("Building {}".format(self.name))
-
+    def build_fundamental_and_external(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -409,7 +406,25 @@ class BaseBatteryModel(pybamm.BaseModel):
             )
             self.variables.update(submodel.get_fundamental_variables())
 
-        # Get coupled variables
+        # set the submodels that are external
+        for sub in self.options["external submodels"]:
+            self.submodels[sub].external = True
+
+        # Set any external variables
+        self.external_variables = []
+        for submodel_name, submodel in self.submodels.items():
+            pybamm.logger.debug(
+                "Getting external variables for {} submodel ({})".format(
+                    submodel_name, self.name
+                )
+            )
+            external_variables = submodel.get_external_variables()
+
+            self.external_variables += external_variables
+
+        self._built_fundamental_and_external = True
+
+    def build_coupled_variables(self):
         # Note: pybamm will try to get the coupled variables for the submodels in the
         # order they are set by the user. If this fails for a particular submodel,
         # return to it later and try again. If setting coupled variables fails and
@@ -444,35 +459,56 @@ class BaseBatteryModel(pybamm.BaseModel):
                             # try setting coupled variables on next loop through
                             pass
 
+    def build_model_equations(self):
         # Set model equations
         for submodel_name, submodel in self.submodels.items():
-            pybamm.logger.debug(
-                "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
-            )
-            submodel.set_rhs(self.variables)
-            pybamm.logger.debug(
-                "Setting algebraic for {} submodel ({})".format(
-                    submodel_name, self.name
+            if submodel.external is False:
+                pybamm.logger.debug(
+                    "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
                 )
-            )
-            submodel.set_algebraic(self.variables)
-            pybamm.logger.debug(
-                "Setting boundary conditions for {} submodel ({})".format(
-                    submodel_name, self.name
+
+                submodel.set_rhs(self.variables)
+                pybamm.logger.debug(
+                    "Setting algebraic for {} submodel ({})".format(
+                        submodel_name, self.name
+                    )
                 )
-            )
-            submodel.set_boundary_conditions(self.variables)
-            pybamm.logger.debug(
-                "Setting initial conditions for {} submodel ({})".format(
-                    submodel_name, self.name
+                submodel.set_algebraic(self.variables)
+                pybamm.logger.debug(
+                    "Setting boundary conditions for {} submodel ({})".format(
+                        submodel_name, self.name
+                    )
                 )
+                submodel.set_boundary_conditions(self.variables)
+                pybamm.logger.debug(
+                    "Setting initial conditions for {} submodel ({})".format(
+                        submodel_name, self.name
+                    )
+                )
+                submodel.set_initial_conditions(self.variables)
+                submodel.set_events(self.variables)
+                pybamm.logger.debug(
+                    "Updating {} submodel ({})".format(submodel_name, self.name)
+                )
+                self.update(submodel)
+
+    def build_model(self):
+
+        # Check if already built
+        if self._built:
+            raise pybamm.ModelError(
+                """Model already built. If you are adding a new submodel, try using
+                `model.update` instead."""
             )
-            submodel.set_initial_conditions(self.variables)
-            submodel.set_events(self.variables)
-            pybamm.logger.debug(
-                "Updating {} submodel ({})".format(submodel_name, self.name)
-            )
-            self.update(submodel)
+
+        pybamm.logger.info("Building {}".format(self.name))
+
+        if self._built_fundamental_and_external is False:
+            self.build_fundamental_and_external()
+
+        self.build_coupled_variables()
+
+        self.build_model_equations()
 
         pybamm.logger.debug("Setting voltage variables")
         self.set_voltage_variables()
@@ -690,6 +726,11 @@ class BaseBatteryModel(pybamm.BaseModel):
             V = pybamm.BoundaryValue(phi_s_cp, "positive tab")
             V_dim = pybamm.BoundaryValue(phi_s_cp_dim, "positive tab")
 
+        phi_s_cn = self.variables["Negative current collector potential"]
+        phi_s_cn_dim = self.variables["Negative current collector potential [V]"]
+        V_local = phi_s_cp - phi_s_cn
+        V_local_dim = phi_s_cp_dim - phi_s_cn_dim
+
         # TODO: add current collector losses to the voltage in 3D
 
         self.variables.update(
@@ -702,6 +743,8 @@ class BaseBatteryModel(pybamm.BaseModel):
                 "X-averaged reaction overpotential [V]": eta_r_av_dim,
                 "X-averaged solid phase ohmic losses": delta_phi_s_av,
                 "X-averaged solid phase ohmic losses [V]": delta_phi_s_av_dim,
+                "Local voltage": V_local,
+                "Local voltage [V]": V_local_dim,
                 "Terminal voltage": V,
                 "Terminal voltage [V]": V_dim,
             }
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index e91d5d676..47291a412 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -46,7 +46,7 @@ class BaseSubModel:
         symbols.
     """
 
-    def __init__(self, param, domain=None, reactions=None):
+    def __init__(self, param, domain=None, reactions=None, external=False):
         super().__init__()
         self.param = param
         # Initialise empty variables (to avoid overwriting with 'None')
@@ -62,6 +62,8 @@ class BaseSubModel:
         self.set_domain_for_broadcast()
         self.reactions = reactions
 
+        self.external = external
+
     @property
     def domain(self):
         return self._domain
@@ -113,6 +115,52 @@ class BaseSubModel:
         """
         return {}
 
+    def get_external_variables(self):
+        """
+        A public method that returns the variables in a submodel which are
+        supplied by an external source.
+
+        Returns
+        -------
+        list :
+            A list of the external variables in the model.
+        """
+
+        external_variables = []
+        list_of_vars = []
+
+        if self.external is True:
+            # look through all the variables in the submodel and get the
+            # variables which are state vectors
+            submodel_variables = self.get_fundamental_variables()
+            for var in submodel_variables.values():
+                if isinstance(var, pybamm.Variable):
+                    list_of_vars += [var]
+
+                elif isinstance(var, pybamm.Concatenation):
+                    if all(
+                        isinstance(child, pybamm.Variable) for child in var.children
+                    ):
+                        list_of_vars += [var]
+
+            # first add only unique concatenations
+            unique_ids = []
+            for var in list_of_vars:
+                if var.id not in unique_ids and isinstance(var, pybamm.Concatenation):
+                    external_variables += [var]
+                    unique_ids += [var.id]
+                    # also add the ids of the children to unique ids
+                    for child in var.children:
+                        unique_ids += [child.id]
+
+            # now add any unique variables that are not part of a concatentation
+            for var in list_of_vars:
+                if var.id not in unique_ids:
+                    external_variables += [var]
+                    unique_ids += [var.id]
+
+        return external_variables
+
     def get_coupled_variables(self, variables):
         """
         A public method that creates and returns the variables in a submodel which
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 1a75fc563..41c383166 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -27,6 +27,26 @@ class Simulation:
     ----------
     model : :class:`pybamm.BaseModel`
         The model to be simulated
+    geometry: :class:`pybamm.Geometry` (optional)
+        The geometry upon which to solve the model
+    parameter_values: dict (optional)
+        A dictionary of parameters and their corresponding numerical
+        values
+    submesh_types: dict (optional)
+        A dictionary of the types of submesh to use on each subdomain
+    var_pts: dict (optional)
+        A dictionary of the number of points used by each spatial
+        variable
+    spatial_methods: dict (optional)
+        A dictionary of the types of spatial method to use on each
+        domain (e.g. pybamm.FiniteVolume)
+    solver: :class:`pybamm.BaseSolver` (optional)
+        The solver to use to solve the model.
+    quick_plot_vars: list (optional)
+        A list of variables to plot automatically
+    C_rate: float (optional)
+        The C_rate at which you would like to run a constant current
+        experiment at.
     """
 
     def __init__(
@@ -39,6 +59,7 @@ class Simulation:
         spatial_methods=None,
         solver=None,
         quick_plot_vars=None,
+        C_rate=None,
     ):
         self.model = model
 
@@ -50,7 +71,13 @@ class Simulation:
         self._solver = solver or self._model.default_solver
         self._quick_plot_vars = quick_plot_vars
 
-        self.reset()
+        self.C_rate = C_rate
+        if self.C_rate:
+            self._parameter_values.update({"C-rate": self.C_rate})
+
+        self._made_first_step = False
+
+        self.reset(update_model=False)
 
         # ignore runtime warnings in notebooks
         if isnotebook():
@@ -71,17 +98,19 @@ class Simulation:
         self._solver = self._model.default_solver
         self._quick_plot_vars = None
 
-    def reset(self):
+    def reset(self, update_model=True):
         """
         A method to reset a simulation back to its unprocessed state.
         """
-        self.model = self.model.new_copy(self._model_options)
+        if update_model:
+            self.model = self.model.new_copy(self._model_options)
         self.geometry = copy.deepcopy(self._unprocessed_geometry)
         self._model_with_set_params = None
         self._built_model = None
         self._mesh = None
         self._disc = None
         self._solution = None
+        self._made_first_step = False
 
     def set_parameters(self):
         """
@@ -147,6 +176,84 @@ class Simulation:
 
         self._solution = solver.solve(self.built_model, t_eval)
 
+    def step(self, dt, solver=None, external_variables=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.
+
+        Parameters
+        ----------
+        dt : numeric type
+            The timestep over which to step the solution
+        solver : :class:`pybamm.BaseSolver`
+            The solver to use to solve the model.
+        external_variables : dict
+            A dictionary of external variables and their corresponding
+            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.
+        save : bool
+            Turn on to store the solution of all previous timesteps
+        """
+        self.build()
+
+        if solver is None:
+            solver = self.solver
+
+        solution = solver.step(
+            self.built_model, dt, external_variables=external_variables
+        )
+
+        if save is False or self._made_first_step is False:
+            self._solution = solution
+        elif self._solution.t[-1] == solution.t[-1]:
+            pass
+        else:
+            self._update_solution(solution)
+
+        self._made_first_step = True
+
+    def _update_solution(self, solution):
+
+        self._solution.set_up_time += solution.set_up_time
+        self._solution.solve_time += solution.solve_time
+        self._solution.t = np.append(self._solution.t, solution.t[-1])
+        self._solution.t_event = solution.t_event
+        self._solution.termination = solution.termination
+        self._solution.y = np.concatenate(
+            [self._solution.y, solution.y[:, -1][:, np.newaxis]], axis=1
+        )
+        self._solution.y_event = solution.y_event
+
+    def get_variable_array(self, *variables):
+        """
+        A helper function to easily obtain a dictionary of arrays of values
+        for a list of variables at the latest timestep.
+
+        Parameters
+        ----------
+        variable: str
+            The name of the variable/variables you wish to obtain the arrays for.
+
+        Returns
+        -------
+        variable_arrays: dict
+            A dictionary of the variable names and their corresponding
+            arrays.
+        """
+
+        variable_arrays = [
+            self.built_model.variables[var].evaluate(
+                self.solution.t[-1], self.solution.y[:, -1]
+            )
+            for var in variables
+        ]
+
+        if len(variable_arrays) == 1:
+            return variable_arrays[0]
+        else:
+            return tuple(variable_arrays)
+
     def plot(self, quick_plot_vars=None, testing=False):
         """
         A method to quickly plot the outputs of the simulation.
@@ -269,6 +376,7 @@ class Simulation:
         spatial_methods=None,
         solver=None,
         quick_plot_vars=None,
+        C_rate=None,
     ):
         """
         A method to set the various specs of the simulation. This method
@@ -291,10 +399,13 @@ class Simulation:
         spatial_methods: dict, optional
             A dictionary of the types of spatial method to use on each
             domain (e.g. pybamm.FiniteVolume)
-        solver: :class:`pybamm.BaseSolver`
+        solver: :class:`pybamm.BaseSolver` (optional)
             The solver to use to solve the model.
-        quick_plot_vars: list
+        quick_plot_vars: list (optional)
             A list of variables to plot automatically
+        C_rate: float (optional)
+            The C_rate at which you would like to run a constant current
+            experiment at.
         """
 
         if model_options:
@@ -316,6 +427,10 @@ class Simulation:
         if quick_plot_vars:
             self._quick_plot_vars = quick_plot_vars
 
+        if C_rate:
+            self.C_rate = C_rate
+            self._parameter_values.update({"C-rate": self.C_rate})
+
         if (
             model_options
             or geometry
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 37f2fd819..5ef78b690 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -22,6 +22,9 @@ class BaseSolver(object):
         self._atol = atol
         self.name = "Base solver"
 
+        self.y_pad = None
+        self.y_ext = None
+
     @property
     def method(self):
         return self._method
@@ -97,7 +100,7 @@ class BaseSolver(object):
         )
         return solution
 
-    def step(self, model, dt, npts=2, log=True):
+    def step(self, model, dt, npts=2, log=True, external_variables=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
@@ -113,6 +116,9 @@ class BaseSolver(object):
         npts : int, optional
             The number of points at which the solution will be returned during
             the step dt. default is 2 (returns the solution at t0 and t0 + dt).
+        external_variables : dict
+            A dictionary of external variables and their corresponding
+            values at the current time
 
         Raises
         ------
@@ -127,6 +133,12 @@ class BaseSolver(object):
         # Set timer
         timer = pybamm.Timer()
 
+        if not hasattr(self, "y0"):
+            # create a y_pad vector of the correct size:
+            self.y_pad = np.zeros((model.y_length - model.external_start, 1))
+
+        self.set_external_variables(model, external_variables)
+
         # Run set up on first step
         if not hasattr(self, "y0"):
             pybamm.logger.info(
@@ -144,22 +156,30 @@ class BaseSolver(object):
                 self.set_up(model)
             self.t = 0.0
             set_up_time = timer.time()
+
         else:
-            set_up_time = None
+            set_up_time = 0
 
         # Step
         t_eval = np.linspace(self.t, self.t + dt, npts)
         solution, solve_time, termination = self.compute_solution(model, t_eval)
 
-        # Assign times
-        solution.solve_time = solve_time
-        if set_up_time:
-            solution.set_up_time = set_up_time
-
         # Set self.t and self.y0 to their values at the final step
         self.t = solution.t[-1]
         self.y0 = solution.y[:, -1]
 
+        # add the external points onto the solution
+        full_y = np.zeros((model.y_length, solution.y.shape[1]))
+        for i in np.arange(solution.y.shape[1]):
+            sol_y = solution.y[:, i]
+            sol_y = sol_y[:, np.newaxis]
+            full_y[:, i] = add_external(sol_y, self.y_pad, self.y_ext)[:, 0]
+        solution.y = full_y
+
+        # Assign times
+        solution.solve_time = solve_time
+        solution.set_up_time = set_up_time
+
         pybamm.logger.debug("Finish stepping {} ({})".format(model.name, termination))
         if set_up_time:
             pybamm.logger.debug(
@@ -175,6 +195,32 @@ class BaseSolver(object):
             )
         return solution
 
+    def set_external_variables(self, model, external_variables):
+        if external_variables is None:
+            external_variables = {}
+
+        # load external variables into a state vector
+        self.y_ext = np.zeros((model.y_length, 1))
+        for var_name, var_vals in external_variables.items():
+            var = model.variables[var_name]
+            if isinstance(var, pybamm.Concatenation):
+                start = var.children[0].y_slices[0].start
+                stop = var.children[-1].y_slices[-1].stop
+                y_slice = slice(start, stop)
+
+            elif isinstance(var, pybamm.StateVector):
+                start = var.y_slices[0].start
+                stop = var.y_slices[-1].stop
+                y_slice = slice(start, stop)
+            else:
+                raise pybamm.InputError(
+                    """The variable you have inputted is not a StateVector or Concatenation
+            of StateVectors. Please check the submodel you have made "external" and
+            ensure that the variable you
+            are passing in is the variable that is solved for in that submodel"""
+                )
+            self.y_ext[y_slice] = var_vals
+
     def compute_solution(self, model, t_eval):
         """Calculate the solution of the model at specified times. Note: this
         does *not* execute the solver setup.
@@ -235,10 +281,21 @@ class BaseSolver(object):
             # Get final event value
             final_event_values = {}
             for name, event in events.items():
+                y_event = add_external(solution.y_event, self.y_pad, self.y_ext)
                 final_event_values[name] = abs(
-                    event.evaluate(solution.t_event, solution.y_event)
+                    event.evaluate(solution.t_event, y_event)
                 )
             termination_event = min(final_event_values, key=final_event_values.get)
             # Add the event to the solution object
             solution.termination = "event: {}".format(termination_event)
             return "the termination event '{}' occurred".format(termination_event)
+
+
+def add_external(y, y_pad, y_ext):
+    """
+    Pad the state vector and then add the external variables so that
+    it is of the correct shape for evaluate
+    """
+    if y_pad is not None and y_ext is not None:
+        y = np.concatenate([y, y_pad]) + y_ext
+    return y
diff --git a/pybamm/solvers/dae_solver.py b/pybamm/solvers/dae_solver.py
index 7ea8283b7..33d312de4 100644
--- a/pybamm/solvers/dae_solver.py
+++ b/pybamm/solvers/dae_solver.py
@@ -7,6 +7,8 @@ import numpy as np
 from scipy import optimize
 from scipy.sparse import issparse
 
+from .base_solver import add_external
+
 
 class DaeSolver(pybamm.BaseSolver):
     """Solve a discretised model.
@@ -79,6 +81,15 @@ class DaeSolver(pybamm.BaseSolver):
         """
         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)
+
         solve_start_time = timer.time()
         pybamm.logger.info("Calling DAE solver")
         solution = self.integrate(
@@ -90,6 +101,7 @@ class DaeSolver(pybamm.BaseSolver):
             jacobian=self.jacobian,
             model=model,
         )
+
         solve_time = timer.time() - solve_start_time
 
         # Identify the event that caused termination
@@ -149,6 +161,8 @@ class DaeSolver(pybamm.BaseSolver):
             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]
 
         else:
@@ -207,58 +221,70 @@ class DaeSolver(pybamm.BaseSolver):
         # 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_alg = casadi.MX.sym(
             "y_alg", len(model.concatenated_algebraic.evaluate(0, y0))
         )
         y_casadi = casadi.vertcat(y_diff, y_alg)
+        if self.y_pad is not None:
+            y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
+            y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+        else:
+            y_casadi_w_ext = y_casadi
+
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
         pybamm.logger.info("Converting algebraic to CasADi")
         concatenated_algebraic = model.concatenated_algebraic.to_casadi(
-            t_casadi, y_casadi
+            t_casadi, y_casadi_w_ext
         )
         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)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext)
             for name, event in model.events.items()
         }
 
         # Create functions to evaluate rhs and algebraic
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
         )
         concatenated_algebraic_fn = casadi.Function(
-            "algebraic", [t_casadi, y_casadi], [concatenated_algebraic]
+            "algebraic", [t_casadi, y_casadi_w_ext], [concatenated_algebraic]
         )
-        all_states_fn = casadi.Function("all", [t_casadi, y_casadi], [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], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
             )
             casadi_jac_alg = casadi.jacobian(concatenated_algebraic, y_casadi)
             casadi_jac_alg_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi], [casadi_jac_alg]
+                "jacobian", [t_casadi, y_casadi_w_ext], [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)
 
         else:
             jacobian = None
             jacobian_alg = None
 
-        # Calculate consistent initial conditions for the algebraic equations
         rhs = RhsCasadi(concatenated_rhs_fn)
         algebraic = AlgebraicCasadi(concatenated_algebraic_fn)
 
+        rhs.set_pad_ext(self.y_pad, self.y_ext)
+        algebraic.set_pad_ext(self.y_pad, self.y_ext)
+
         if len(model.algebraic) > 0:
 
             y0 = self.calculate_consistent_initial_conditions(
@@ -273,7 +299,9 @@ class DaeSolver(pybamm.BaseSolver):
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
+            casadi_event_fn = casadi.Function(
+                "event", [t_casadi, y_casadi_w_ext], [event]
+            )
             return EvalEvent(casadi_event_fn)
 
         # Add the solver attributes
@@ -416,8 +444,16 @@ class Rhs:
 
     def __init__(self, concatenated_rhs_fn):
         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):
+        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]
 
 
@@ -425,6 +461,8 @@ class RhsCasadi(Rhs):
     "Returns information about rhs 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.concatenated_rhs_fn(t, y).full()[:, 0]
 
 
@@ -433,8 +471,16 @@ class Algebraic:
 
     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]
 
 
@@ -442,6 +488,8 @@ class AlgebraicCasadi(Algebraic):
     "Returns information about algebraic equations 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.concatenated_algebraic_fn(t, y).full()[:, 0]
 
 
@@ -453,12 +501,19 @@ class Residuals:
         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(
             "Evaluating residuals for {} at t={}".format(self.model.name, t)
         )
         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={})
         # reuse known_evals
         alg_eval = self.concatenated_algebraic_fn(t, y, known_evals=known_evals)[0]
@@ -475,11 +530,15 @@ class ResidualsCasadi(Residuals):
         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(
             "Evaluating residuals for {} at t={}".format(self.model.name, t)
         )
+        y = y[:, np.newaxis]
+        y = add_external(y, self.y_pad, self.y_ext)
         states_eval = self.all_states_fn(t, y).full()[:, 0]
         return states_eval - self.mass_matrix @ ydot
 
@@ -489,8 +548,16 @@ class EvalEvent:
 
     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)
 
 
@@ -499,8 +566,16 @@ class Jacobian:
 
     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]
 
 
@@ -508,4 +583,7 @@ class JacobianCasadi(Jacobian):
     "Returns information about 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)
+
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index e193341c8..6c16cf008 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -260,8 +260,13 @@ class IDAKLUSolver(pybamm.DaeSolver):
             # 2 = found root(s)
             elif sol.flag == 2:
                 termination = "event"
+
             return pybamm.Solution(
-                sol.t, np.transpose(y_out), t[-1], np.transpose(y_out[-1]), termination
+                sol.t,
+                np.transpose(y_out),
+                t[-1],
+                np.transpose(y_out[-1])[:, np.newaxis],
+                termination,
             )
         else:
             raise pybamm.SolverError(sol.message)
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 6feba0d98..934e521ef 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -5,6 +5,8 @@ import casadi
 import pybamm
 import numpy as np
 
+from .base_solver import add_external
+
 
 class OdeSolver(pybamm.BaseSolver):
     """Solve a discretised model.
@@ -35,6 +37,12 @@ class OdeSolver(pybamm.BaseSolver):
         """
         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)
+
         solve_start_time = timer.time()
         pybamm.logger.info("Calling ODE solver")
         solution = self.integrate(
@@ -45,6 +53,7 @@ class OdeSolver(pybamm.BaseSolver):
             mass_matrix=model.mass_matrix.entries,
             jacobian=self.jacobian,
         )
+
         solve_time = timer.time() - solve_start_time
 
         # Identify the event that caused termination
@@ -163,22 +172,31 @@ class OdeSolver(pybamm.BaseSolver):
 
         t_casadi = casadi.MX.sym("t")
         y_casadi = casadi.MX.sym("y", len(y0))
+
+        if self.y_pad is not None:
+            y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
+            y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
+        else:
+            y_casadi_w_ext = y_casadi
+
         pybamm.logger.info("Converting RHS to CasADi")
-        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi)
+        concatenated_rhs = model.concatenated_rhs.to_casadi(t_casadi, y_casadi_w_ext)
         pybamm.logger.info("Converting events to CasADi")
         casadi_events = {
-            name: event.to_casadi(t_casadi, y_casadi)
+            name: event.to_casadi(t_casadi, y_casadi_w_ext)
             for name, event in model.events.items()
         }
 
         # Create function to evaluate rhs
         concatenated_rhs_fn = casadi.Function(
-            "rhs", [t_casadi, y_casadi], [concatenated_rhs]
+            "rhs", [t_casadi, y_casadi_w_ext], [concatenated_rhs]
         )
 
         # Create event-dependent function to evaluate events
         def get_event_class(event):
-            casadi_event_fn = casadi.Function("event", [t_casadi, y_casadi], [event])
+            casadi_event_fn = casadi.Function(
+                "event", [t_casadi, y_casadi_w_ext], [event]
+            )
             return EvalEvent(casadi_event_fn)
 
         # Create function to evaluate jacobian
@@ -186,8 +204,9 @@ class OdeSolver(pybamm.BaseSolver):
             pybamm.logger.info("Calculating jacobian")
             casadi_jac = casadi.jacobian(concatenated_rhs, y_casadi)
             casadi_jac_fn = casadi.Function(
-                "jacobian", [t_casadi, y_casadi], [casadi_jac]
+                "jacobian", [t_casadi, y_casadi_w_ext], [casadi_jac]
             )
+
             jacobian = JacobianCasadi(casadi_jac_fn)
 
         else:
@@ -226,18 +245,23 @@ class OdeSolver(pybamm.BaseSolver):
 
 
 # Set up caller classes outside of the solver object to allow pickling
-
-
 class Dydt:
     "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]
         return dy[:, 0]
 
@@ -247,6 +271,8 @@ class DydtCasadi(Dydt):
 
     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()
         return dy[:, 0]
 
@@ -256,8 +282,16 @@ class EvalEvent:
 
     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)
 
 
@@ -266,8 +300,16 @@ class Jacobian:
 
     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]
 
 
@@ -275,4 +317,6 @@ class JacobianCasadi(Jacobian):
     "Returns information about 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)
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index e0e8f3970..98f8898f3 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -86,6 +86,32 @@ class FiniteVolume(pybamm.SpatialMethod):
         out = gradient_matrix @ discretised_symbol
         return out
 
+    def preprocess_external_variables(self, var):
+        """
+        For finite volumes, we need the boundary fluxes for discretising
+        properly. Here, we extrapolate and then add them to the boundary
+        conditions.
+
+        Parameters
+        ----------
+        var : :class:`pybamm.Variable` or :class:`pybamm.Concatenation`
+            The external variable that is to be processed
+
+        Returns
+        -------
+        new_bcs: dict
+            A dictionary containing the new boundary conditions
+        """
+
+        new_bcs = {
+            var: {
+                "left": (pybamm.BoundaryGradient(var, "left"), "Neumann"),
+                "right": (pybamm.BoundaryGradient(var, "right"), "Neumann"),
+            }
+        }
+
+        return new_bcs
+
     def gradient_matrix(self, domain):
         """
         Gradient matrix for finite volumes in the appropriate domain.
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 0ed9e70dc..1258bc0f9 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -292,6 +292,9 @@ class SpatialMethod:
 
         raise NotImplementedError
 
+    def preprocess_external_variables(self, var):
+        return {}
+
     def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         """
         Returns the boundary value or flux using the approriate expression for the
diff --git a/results/2plus1D/set_temperature_spm_1plus1D.py b/results/2plus1D/set_temperature_spm_1plus1D.py
index f6101a882..c469d7607 100644
--- a/results/2plus1D/set_temperature_spm_1plus1D.py
+++ b/results/2plus1D/set_temperature_spm_1plus1D.py
@@ -4,189 +4,48 @@
 
 import pybamm
 import numpy as np
-import matplotlib.pyplot as plt
-import sys
 
 # set logging level
 pybamm.set_logging_level("INFO")
 
-# load (1+1D) SPM model
-options = {
+model_options = {
     "current collector": "potential pair",
     "dimensionality": 1,
-    "thermal": "set external temperature",
+    "thermal": "x-lumped",
+    "external submodels": ["thermal"],
 }
-model = pybamm.lithium_ion.SPM(options)
+model = pybamm.lithium_ion.SPMe(model_options)
 
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
-
-# set mesh
-nbat = 5
 var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: nbat}
-# depending on number of points in y-z plane may need to increase recursion depth...
-sys.setrecursionlimit(10000)
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
-
-
-# define a method which updates statevector
-def update_statevector(variables, statevector):
-    "takes in a dict of variable name and vector of updated state"
-    for name, new_vector in variables.items():
-        var_slice = model.variables[name].y_slices
-        statevector[var_slice] = new_vector
-    return statevector
-
-
-# define a method to non-dimensionalise a temperature
-def non_dim_temperature(temperature):
-    "takes in a temperature and returns the non-dimensional version"
-    Delta_T = param.process_symbol(model.submodels["thermal"].param.Delta_T).evaluate()
-    T_ref = param.process_symbol(model.submodels["thermal"].param.T_ref).evaluate()
-    return (temperature - T_ref) / Delta_T
-
-
-# step model in time and process variables for later plotting
-solver = model.default_solver
-dt = 0.1  # timestep to take
-npts = 20  # number of points to store the solution at during this step
-solution1 = solver.step(model, dt, npts=npts)
-# create dict of variables to post process
-output_variables = [
-    "Negative current collector potential [V]",
-    "Positive current collector potential [V]",
-    "Current [A]",
-    "X-averaged total heating [W.m-3]",
-    "X-averaged positive particle surface concentration [mol.m-3]",
-    "X-averaged cell temperature [K]",
-]
-output_variables_dict = {}
-for var in output_variables:
-    output_variables_dict[var] = model.variables[var]
-processed_vars_step1 = pybamm.post_process_variables(
-    output_variables_dict, solution1.t, solution1.y, mesh
+z_pts = 20
+var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 5, var.r_p: 5, var.z: z_pts}
+
+sim = pybamm.Simulation(model, var_pts=var_pts, C_rate=2)
+
+# Set the temperature (in dimensionless form)
+# T_av = np.linspace(0, 1, z_pts)[:, np.newaxis]
+
+z = np.linspace(0, 1, z_pts)
+t_eval = np.linspace(0, 0.13, 50)
+# step through the solver, setting the temperature at each timestep
+for i in np.arange(1, len(t_eval) - 1):
+    dt = t_eval[i + 1] - t_eval[i]
+    T_av = (np.sin(2 * np.pi * z) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+        :, np.newaxis
+    ]
+    external_variables = {"X-averaged cell temperature": T_av}
+    sim.step(dt, external_variables=external_variables)
+
+sim.plot(
+    [
+        "Terminal voltage [V]",
+        "X-averaged total heating [W.m-3]",
+        "X-averaged cell temperature [K]",
+        "X-averaged negative particle surface concentration [mol.m-3]",
+        "X-averaged positive particle surface concentration [mol.m-3]",
+        "Negative current collector potential [V]",
+        "Positive current collector potential [V]",
+        "Local voltage [V]",
+    ]
 )
 
-# get the current state and temperature
-current_state = solution1.y[:, -1]
-temp_ave = current_state[model.variables["X-averaged cell temperature"].y_slices]
-
-# update the temperature
-T_ref = param.process_symbol(model.submodels["thermal"].param.T_ref).evaluate()
-t_external = np.linspace(T_ref, T_ref + 6.0, nbat)
-non_dim_t_external = non_dim_temperature(t_external)
-variables = {"X-averaged cell temperature": non_dim_t_external}
-new_state = update_statevector(variables, current_state)
-
-# step in time again and process variables for later plotting
-# use new state as initial condition. Note: need to to recompute consistent initial
-# values for the algebraic part of the model. Since the (dummy) equation for the
-# temperature is an ODE, the imposed change in temperature is unaffected by this
-# process
-solver.y0 = solver.calculate_consistent_initial_conditions(
-    solver.rhs, solver.algebraic, new_state
-)
-solution2 = solver.step(model, dt, npts=npts)
-processed_vars_step2 = pybamm.post_process_variables(
-    output_variables_dict, solution2.t, solution2.y, mesh
-)
-
-# plots
-t_sec = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge)
-t_hour = t_sec / (3600)
-z = np.linspace(0, 1, nbat)
-
-# local voltage
-plt.figure()
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1["Positive current collector potential [V]"](
-            solution1.t, z=z
-        )[bat_id, :]
-        - processed_vars_step1["Negative current collector potential [V]"](
-            solution1.t, z=z
-        )[bat_id, :],
-        solution2.t * t_hour,
-        processed_vars_step2["Positive current collector potential [V]"](
-            solution2.t, z=z
-        )[bat_id, :]
-        - processed_vars_step1["Negative current collector potential [V]"](
-            solution2.t, z=z
-        )[bat_id, :],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("Local voltage [V]")
-
-# applied current
-plt.figure()
-plt.plot(
-    solution1.t,
-    processed_vars_step1["Current [A]"](solution1.t),
-    solution2.t,
-    processed_vars_step2["Current [A]"](solution2.t),
-)
-plt.xlabel("t")
-plt.ylabel("Current [A]")
-
-# local heating
-plt.figure()
-z = np.linspace(0, 1, nbat)
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1["X-averaged total heating [W.m-3]"](solution1.t, z=z)[
-            bat_id, :
-        ],
-        solution2.t * t_hour,
-        processed_vars_step2["X-averaged total heating [W.m-3]"](solution2.t, z=z)[
-            bat_id, :
-        ],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("X-averaged total heating [W.m-3]")
-plt.yscale("log")
-
-# local concentration
-plt.figure()
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1[
-            "X-averaged positive particle surface concentration [mol.m-3]"
-        ](solution1.t, z=z)[bat_id, :],
-        solution2.t * t_hour,
-        processed_vars_step2[
-            "X-averaged positive particle surface concentration [mol.m-3]"
-        ](solution2.t, z=z)[bat_id, :],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("X-averaged positive particle surface concentration [mol.m-3]")
-
-# local temperature
-plt.figure()
-for bat_id in range(nbat):
-    plt.plot(
-        solution1.t * t_hour,
-        processed_vars_step1["X-averaged cell temperature [K]"](solution1.t, z=z)[
-            bat_id, :
-        ],
-        solution2.t * t_hour,
-        processed_vars_step2["X-averaged cell temperature [K]"](solution2.t, z=z)[
-            bat_id, :
-        ],
-    )
-plt.xlabel("t [hrs]")
-plt.ylabel("X-averaged cell temperature [K]")
-
-plt.show()


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/__init__.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
new file mode 100644
index 000000000..7f0a661fd
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
@@ -0,0 +1,58 @@
+#
+# Tests for external submodels
+#
+import pybamm
+import unittest
+import numpy as np
+
+
+class TestExternalCC(unittest.TestCase):
+    @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
+    def test_2p1d(self):
+        model_options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "external submodels": ["current collector"],
+        }
+        model = pybamm.lithium_ion.DFN(model_options)
+        yz_pts = 3
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: 4,
+            pybamm.standard_spatial_vars.x_s: 4,
+            pybamm.standard_spatial_vars.x_p: 4,
+            pybamm.standard_spatial_vars.r_n: 4,
+            pybamm.standard_spatial_vars.r_p: 4,
+            pybamm.standard_spatial_vars.y: yz_pts,
+            pybamm.standard_spatial_vars.z: yz_pts,
+        }
+        solver = pybamm.IDAKLUSolver()
+        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
+
+        t_eval = np.linspace(0, 0.08, 3)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+
+            # provide phi_s_n and i_cc
+            phi_s_n = np.zeros((yz_pts ** 2, 1))
+            i_boundary_cc = np.ones((yz_pts ** 2, 1))
+            external_variables = {
+                "Negative current collector potential": phi_s_n,
+                "Current collector current density": i_boundary_cc,
+            }
+
+            sim.step(dt, external_variables=external_variables)
+
+            # obtain phi_s_n from the pybamm solution at the current time
+            phi_s_p = sim.get_variable_array("Positive current collector potential")
+
+        self.assertTrue(phi_s_p.shape, (yz_pts ** 2, 1))
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
new file mode 100644
index 000000000..26600fc46
--- /dev/null
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_temperature.py
@@ -0,0 +1,109 @@
+#
+# Tests for external submodels
+#
+import pybamm
+import unittest
+import numpy as np
+
+
+class TestExternalThermalModels(unittest.TestCase):
+    def test_external_lumped_temperature(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+
+        t_eval = np.linspace(0, 0.01, 3)
+
+        T_av = 0
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            external_variables = {"X-averaged cell temperature": T_av}
+            T_av += 1
+            sim.step(dt, external_variables=external_variables)
+
+    @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
+    def test_external_temperature(self):
+
+        model_options = {
+            "thermal": "x-full",
+            "external submodels": ["thermal"],
+        }
+
+        model = pybamm.lithium_ion.SPM(model_options)
+
+        neg_pts = 5
+        sep_pts = 3
+        pos_pts = 5
+        tot_pts = neg_pts + sep_pts + pos_pts
+
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: neg_pts,
+            pybamm.standard_spatial_vars.x_s: sep_pts,
+            pybamm.standard_spatial_vars.x_p: pos_pts,
+            pybamm.standard_spatial_vars.r_n: 5,
+            pybamm.standard_spatial_vars.r_p: 5,
+        }
+
+        sim = pybamm.Simulation(model, var_pts=var_pts)
+
+        t_eval = np.linspace(0, 0.01, 3)
+        x = np.linspace(0, 1, tot_pts)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+                :, np.newaxis
+            ]
+            external_variables = {"Cell temperature": T}
+            sim.step(dt, external_variables=external_variables)
+
+    @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
+    def test_dae_external_temperature(self):
+
+        model_options = {
+            "thermal": "x-full",
+            "external submodels": ["thermal"],
+        }
+
+        model = pybamm.lithium_ion.DFN(model_options)
+
+        neg_pts = 5
+        sep_pts = 3
+        pos_pts = 5
+        tot_pts = neg_pts + sep_pts + pos_pts
+
+        var_pts = {
+            pybamm.standard_spatial_vars.x_n: neg_pts,
+            pybamm.standard_spatial_vars.x_s: sep_pts,
+            pybamm.standard_spatial_vars.x_p: pos_pts,
+            pybamm.standard_spatial_vars.r_n: 5,
+            pybamm.standard_spatial_vars.r_p: 5,
+        }
+
+        solver = pybamm.IDAKLUSolver()
+        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
+        sim.build()
+
+        t_eval = np.linspace(0, 0.01, 3)
+        x = np.linspace(0, 1, tot_pts)
+
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            T = (np.sin(2 * np.pi * x) * np.sin(2 * np.pi * 100 * t_eval[i]))[
+                :, np.newaxis
+            ]
+            external_variables = {"Cell temperature": T}
+            sim.step(dt, external_variables=external_variables)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index 7e56b605e..ec064b75c 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -16,6 +16,7 @@ class TestSPM(unittest.TestCase):
 
     def test_basic_processing_1plus1D(self):
         options = {"current collector": "potential pair", "dimensionality": 1}
+
         model = pybamm.lithium_ion.SPM(options)
         var = pybamm.standard_spatial_vars
         var_pts = {
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 0a61b0507..b664f1e43 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -66,6 +66,68 @@ class TestDiscretise(unittest.TestCase):
         for child in c_e.children:
             self.assertTrue(child.id in disc.bcs.keys())
 
+    def test_adding_0D_external_variable(self):
+        model = pybamm.BaseModel()
+        a = pybamm.Variable("a")
+        b = pybamm.Variable("b")
+
+        model.rhs = {a: a * b}
+        model.initial_conditions = {a: 0}
+        model.external_variables = [b]
+        model.variables = {"a": a, "b": b, "c": a * b}
+
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        self.assertEqual(len(model.y_slices), 2)
+        self.assertEqual(model.y_slices[b.id][0], slice(1, 2, None))
+
+    def test_adding_1D_external_variable(self):
+        model = pybamm.BaseModel()
+
+        a = pybamm.Variable("a", domain=["test"])
+        b = pybamm.Variable("b", domain=["test"])
+
+        model.rhs = {a: a * b}
+        model.boundary_conditions = {
+            a: {"left": (0, "Dirichlet"), "right": (0, "Dirichlet")}
+        }
+        model.initial_conditions = {a: 0}
+        model.external_variables = [b]
+        model.variables = {
+            "a": a,
+            "b": b,
+            "c": a * b,
+            "grad b": pybamm.grad(b),
+            "div grad b": pybamm.div(pybamm.grad(b)),
+        }
+
+        x = pybamm.SpatialVariable("x", domain="test", coord_sys="cartesian")
+        geometry = {
+            "test": {"primary": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
+        }
+
+        submesh_types = {"test": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh)}
+        var_pts = {x: 10}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+        spatial_methods = {"test": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        self.assertEqual(len(model.y_slices), 2)
+        self.assertEqual(model.y_slices[a.id][0], slice(0, 10, None))
+        self.assertEqual(model.y_slices[b.id][0], slice(10, 20, None))
+
+        # check that b is added to the boundary conditions
+        model.bcs[b.id]["left"]
+        model.bcs[b.id]["right"]
+
+        # check that grad and div(grad ) produce the correct shapes
+        self.assertEqual(model.variables["b"].shape_for_testing, (10, 1))
+        self.assertEqual(model.variables["grad b"].shape_for_testing, (11, 1))
+        self.assertEqual(model.variables["div grad b"].shape_for_testing, (10, 1))
+
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()
diff --git a/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py b/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py
index 033fb95d8..f0d45ae51 100644
--- a/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py
+++ b/tests/unit/test_models/test_submodels/standard_submodel_unit_tests.py
@@ -14,10 +14,15 @@ class StandardSubModelTests(object):
             variables = {}
         self.submodel = submodel
         self.variables = variables
+        self.external_variables = []
 
     def test_get_fundamental_variables(self):
         self.variables.update(self.submodel.get_fundamental_variables())
 
+    def test_get_external_variables(self):
+        external_variables = self.submodel.get_external_variables()
+        self.external_variables += external_variables
+
     def test_get_coupled_variables(self):
         self.variables.update(self.submodel.get_coupled_variables(self.variables))
 
@@ -38,6 +43,7 @@ class StandardSubModelTests(object):
 
     def test_all(self):
         self.test_get_fundamental_variables()
+        self.test_get_external_variables()
         self.test_get_coupled_variables()
         self.test_set_rhs()
         self.test_set_algebraic()
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 25cd6c7bf..21755dbd7 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -1,6 +1,6 @@
 import pybamm
-import unittest
 import numpy as np
+import unittest
 
 
 class TestSimulation(unittest.TestCase):
@@ -130,6 +130,89 @@ class TestSimulation(unittest.TestCase):
         sim.specs(spatial_methods=spatial_methods)
         sim.build()
 
+    def test_set_crate(self):
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM(), C_rate=2)
+        self.assertEqual(sim.parameter_values["C-rate"], 2)
+        sim.specs(C_rate=3)
+        self.assertEqual(sim.parameter_values["C-rate"], 3)
+
+    def test_set_defaults(self):
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+
+        model_options = {"thermal": "x-full"}
+        submesh_types = {
+            "Negative particle": pybamm.MeshGenerator(pybamm.Exponential1DSubMesh)
+        }
+        solver = pybamm.BaseSolver()
+        quick_plot_vars = ["Negative particle surface concentration"]
+        sim.specs(
+            model_options=model_options,
+            submesh_types=submesh_types,
+            solver=solver,
+            quick_plot_vars=quick_plot_vars,
+        )
+
+        sim.set_defaults()
+
+        self.assertEqual(sim.model_options["thermal"], "x-full")
+        self.assertEqual(
+            sim.submesh_types["negative particle"].submesh_type, pybamm.Uniform1DSubMesh
+        )
+        self.assertEqual(sim.quick_plot_vars, None)
+        self.assertIsInstance(sim.solver, pybamm.ScipySolver)
+
+    def test_get_variable_array(self):
+
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+        sim.solve()
+
+        phi_s_n = sim.get_variable_array("Negative electrode potential")
+
+        self.assertIsInstance(phi_s_n, np.ndarray)
+
+        c_s_n_surf, c_e = sim.get_variable_array(
+            "Negative particle surface concentration", "Electrolyte concentration"
+        )
+
+        self.assertIsInstance(c_s_n_surf, np.ndarray)
+        self.assertIsInstance(c_e, np.ndarray)
+
+    def test_set_external_variable(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+
+        T_av = 0
+
+        dt = 0.001
+
+        external_variables = {"X-averaged cell temperature": T_av}
+        sim.step(dt, external_variables=external_variables)
+
+    def test_step(self):
+
+        dt = 0.001
+        sim = pybamm.Simulation(pybamm.lithium_ion.SPM())
+        sim.step(dt)  # 1 step stores first two points
+        self.assertEqual(sim.solution.t.size, 2)
+        self.assertEqual(sim.solution.y[0, :].size, 2)
+        self.assertEqual(sim.solution.t[0], 0)
+        self.assertEqual(sim.solution.t[1], dt)
+        sim.step(dt)  # automatically append the next step
+        self.assertEqual(sim.solution.t.size, 3)
+        self.assertEqual(sim.solution.y[0, :].size, 3)
+        self.assertEqual(sim.solution.t[0], 0)
+        self.assertEqual(sim.solution.t[1], dt)
+        self.assertEqual(sim.solution.t[2], 2 * dt)
+        sim.step(dt, save=False)  # now only store the two end step points
+        self.assertEqual(sim.solution.t.size, 2)
+        self.assertEqual(sim.solution.y[0, :].size, 2)
+        self.assertEqual(sim.solution.t[0], 2 * dt)
+        self.assertEqual(sim.solution.t[1], 3 * dt)
+
     def test_save_load(self):
         model = pybamm.lead_acid.LOQS()
         model.use_jacobian = True
@@ -193,7 +276,7 @@ class TestSimulation(unittest.TestCase):
         sim_load = pybamm.load_sim("test.pickle")
         self.assertEqual(sim.model.name, sim_load.model.name)
 
-    def test_set_defaults(self):
+    def test_set_defaults2(self):
         model = pybamm.lithium_ion.SPM()
 
         # make simulation with silly options (should this be allowed?)
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 9d0075e3e..01fb85bcc 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -2,6 +2,7 @@
 # Tests for the Base Solver class
 #
 import pybamm
+import numpy as np
 
 import unittest
 
@@ -30,6 +31,17 @@ class TestBaseSolver(unittest.TestCase):
         with self.assertRaisesRegex(pybamm.ModelError, "Cannot solve empty model"):
             solver.solve(model, None)
 
+    def test_set_external_variables(self):
+        options = {"thermal": "x-full", "external submodels": ["thermal"]}
+        model = pybamm.lithium_ion.SPM(options)
+        sim = pybamm.Simulation(model)
+        sim.build()
+        solver = pybamm.BaseSolver()
+
+        T = np.ones((60, 1))
+        external_variables = {"Cell temperature": T}
+        solver.set_external_variables(sim.built_model, external_variables)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")