diff --git a/CHANGELOG.md b/CHANGELOG.md
index 602ed80b21..11424767f9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@
 
 ## Features
 
-- Added variables "Loss of lithium due to LAM in negative/positive electrode [mol]", and "Total lithium lost to LAM [mol]/[A.h]". This is also included in the calculation of "total lithium in system" to make sure that lithium is truly conserved. ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
+- Added variables "Loss of lithium due to loss of active material in negative/positive electrode [mol]". These should be included in the calculation of "total lithium in system" to make sure that lithium is truly conserved. ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
 - `initial_soc` can now be a string "x V", in which case the simulation is initialized to start from that voltage ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 - The `ElectrodeSOH` solver can now calculate electrode balance based on a target "cell capacity" (requires cell capacity "Q" as input), as well as the default "cyclable cell capacity" (requires cyclable lithium capacity "Q_Li" as input). Use the keyword argument `known_value` to control which is used. ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 
@@ -14,6 +14,7 @@
 ## Breaking changes
 
 - Renamed "Negative/Positive electrode SOC" to "Negative/Positive electrode stoichiometry" to avoid confusion with cell SOC ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
+- Removed external variables and submodels. InputParameter should now be used in all cases ([#2502](https://github.com/pybamm-team/PyBaMM/pull/2502))
 - Trying to use a solver to solve multiple models results in a RuntimeError exception ([#2481](https://github.com/pybamm-team/PyBaMM/pull/2481))
 - Inputs for the `ElectrodeSOH` solver are now (i) "Q_Li", the total cyclable capacity of lithium in the electrodes (previously "n_Li", the total number of moles, n_Li = 3600/F \* Q_Li) (ii) "Q_n", the capacity of the negative electrode (previously "C_n"), and "Q_p", the capacity of the positive electrode (previously "C_p") ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 
diff --git a/docs/source/expression_tree/variable.rst b/docs/source/expression_tree/variable.rst
index 510e366cf1..ab42f76bcf 100644
--- a/docs/source/expression_tree/variable.rst
+++ b/docs/source/expression_tree/variable.rst
@@ -7,6 +7,3 @@ Variable
 .. autoclass:: pybamm.VariableDot
   :members:
 
-.. autoclass:: pybamm.ExternalVariable
-  :members:
-
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index f2dff27cd7..d2e50d616a 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -63,7 +63,6 @@ def __init__(self, mesh=None, spatial_methods=None):
         self._bcs = {}
         self.y_slices = {}
         self._discretised_symbols = {}
-        self.external_variables = {}
 
     @property
     def mesh(self):
@@ -178,11 +177,6 @@ def process_model(
         pybamm.logger.verbose("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)
-        self.set_external_variables(model)
-
         # set boundary conditions (only need key ids for boundary_conditions)
         pybamm.logger.verbose(
             "Discretise boundary conditions for {}".format(model.name)
@@ -242,11 +236,6 @@ def process_model(
             processed_events.append(processed_event)
         model_disc.events = processed_events
 
-        # Set external variables
-        model_disc.external_variables = [
-            self.process_symbol(var) for var in model.external_variables
-        ]
-
         # Create mass matrix
         pybamm.logger.verbose("Create mass matrix for {}".format(model.name))
         model_disc.mass_matrix, model_disc.mass_matrix_inv = self.create_mass_matrix(
@@ -349,70 +338,6 @@ def _get_variable_size(self, variable):
                 size += spatial_method.mesh[dom].npts_for_broadcast_to_nodes * repeats
             return size
 
-    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_external_variables(self, model):
-        """
-        Add external variables to the list of variables to account for, being careful
-        about concatenations
-        """
-        for var in model.external_variables:
-            # Find the name in the model variables
-            # Look up dictionary key based on value
-            try:
-                idx = list(model.variables.values()).index(var)
-            except ValueError:
-                raise ValueError(
-                    """
-                    Variable {} must be in the model.variables dictionary to be set
-                    as an external variable
-                    """.format(
-                        var
-                    )
-                )
-            name = list(model.variables.keys())[idx]
-            if isinstance(var, pybamm.Variable):
-                # No need to keep track of the parent
-                self.external_variables[(name, None)] = var
-            elif isinstance(var, pybamm.ConcatenationVariable):
-                start = 0
-                end = 0
-                for child in var.children:
-                    dom = child.domain[0]
-                    if (
-                        self.spatial_methods[dom]._get_auxiliary_domain_repeats(
-                            child.domains
-                        )
-                        > 1
-                    ):
-                        raise NotImplementedError(
-                            "Cannot create 2D external variable with concatenations"
-                        )
-                    end += self._get_variable_size(child)
-                    # Keep a record of the parent
-                    self.external_variables[(name, (var, start, end))] = child
-                    start = end
-
     def set_internal_boundary_conditions(self, model):
         """
         A method to set the internal boundary conditions for the submodel.
@@ -526,16 +451,15 @@ def process_boundary_conditions(self, model):
 
             # check if the boundary condition at the origin for sphere domains is other
             # than no flux
-            if key not in model.external_variables:
-                for subdomain in key.domain:
-                    if self.mesh[subdomain].coord_sys == "spherical polar":
-                        if bcs["left"][0].value != 0 or bcs["left"][1] != "Neumann":
-                            raise pybamm.ModelError(
-                                "Boundary condition at r = 0 must be a homogeneous "
-                                "Neumann condition for {} coordinates".format(
-                                    self.mesh[subdomain].coord_sys
-                                )
+            for subdomain in key.domain:
+                if self.mesh[subdomain].coord_sys == "spherical polar":
+                    if bcs["left"][0].value != 0 or bcs["left"][1] != "Neumann":
+                        raise pybamm.ModelError(
+                            "Boundary condition at r = 0 must be a homogeneous "
+                            "Neumann condition for {} coordinates".format(
+                                self.mesh[subdomain].coord_sys
                             )
+                        )
 
             # Handle any boundary conditions applied on the tabs
             if any("tab" in side for side in list(bcs.keys())):
@@ -956,52 +880,26 @@ def _process_symbol(self, symbol):
             )
 
         elif isinstance(symbol, pybamm.Variable):
-            # Check if variable is a standard variable or an external variable
-            if any(symbol == var for var in self.external_variables.values()):
-                # Look up dictionary key based on value
-                idx = list(self.external_variables.values()).index(symbol)
-                name, parent_and_slice = list(self.external_variables.keys())[idx]
-                if parent_and_slice is None:
-                    # Variable didn't come from a concatenation so we can just create a
-                    # normal external variable using the symbol's name
-                    return pybamm.ExternalVariable(
-                        symbol.name,
-                        size=self._get_variable_size(symbol),
-                        domains=symbol.domains,
-                    )
-                else:
-                    # We have to use a special name since the concatenation doesn't have
-                    # a very informative name. Needs improving
-                    parent, start, end = parent_and_slice
-                    ext = pybamm.ExternalVariable(
-                        name,
-                        size=self._get_variable_size(parent),
-                        domains=parent.domains,
-                    )
-                    out = pybamm.Index(ext, slice(start, end))
-                    out.copy_domains(symbol)
-                    return out
-            else:
-                # add a try except block for a more informative error if a variable
-                # can't be found. This should usually be caught earlier by
-                # model.check_well_posedness, but won't be if debug_mode is False
-                try:
-                    y_slices = self.y_slices[symbol]
-                except KeyError:
-                    raise pybamm.ModelError(
-                        """
-                        No key set for variable '{}'. Make sure it is included in either
-                        model.rhs, model.algebraic, or model.external_variables in an
-                        unmodified form (e.g. not Broadcasted)
-                        """.format(
-                            symbol.name
-                        )
+            # add a try except block for a more informative error if a variable
+            # can't be found. This should usually be caught earlier by
+            # model.check_well_posedness, but won't be if debug_mode is False
+            try:
+                y_slices = self.y_slices[symbol]
+            except KeyError:
+                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)
+                    """.format(
+                        symbol.name
                     )
-                # Add symbol's reference and multiply by the symbol's scale
-                # so that the state vector is of order 1
-                return symbol.reference + symbol.scale * pybamm.StateVector(
-                    *y_slices, domains=symbol.domains
                 )
+            # Add symbol's reference and multiply by the symbol's scale
+            # so that the state vector is of order 1
+            return symbol.reference + symbol.scale * pybamm.StateVector(
+                *y_slices, domains=symbol.domains
+            )
 
         elif isinstance(symbol, pybamm.SpatialVariable):
             return spatial_method.spatial_variable(symbol)
@@ -1083,14 +981,7 @@ def _concatenate_in_order(self, var_eqn_dict, check_complete=False, sparse=False
         if check_complete:
             # Check keys from the given var_eqn_dict against self.y_slices
             unpacked_variables_set = set(unpacked_variables)
-            external_vars = set(self.external_variables.values())
-            for var in self.external_variables.values():
-                child_vars = set(var.children)
-                external_vars = external_vars.union(child_vars)
-            y_slices_with_external_removed = set(self.y_slices.keys()).difference(
-                external_vars
-            )
-            if unpacked_variables_set != y_slices_with_external_removed:
+            if unpacked_variables_set != set(self.y_slices.keys()):
                 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 1a66d230c3..3d064df520 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -56,12 +56,6 @@ class ModelWarning(UserWarning):
     pass
 
 
-class InputError(Exception):
-    """An external variable has been input incorrectly into PyBaMM."""
-
-    pass
-
-
 class DiscretisationError(Exception):
     """A model could not be discretised."""
 
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index 661771a123..b3a048b1f1 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -55,7 +55,6 @@ def _convert(self, symbol, t, y, y_dot, inputs):
                 pybamm.Array,
                 pybamm.Time,
                 pybamm.InputParameter,
-                pybamm.ExternalVariable,
             ),
         ):
             return casadi.MX(symbol.evaluate(t, y, y_dot, inputs))
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 61f6ed8f55..836bd81bca 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -1,7 +1,6 @@
 #
 # Variable class
 #
-import numbers
 
 import numpy as np
 import sympy
@@ -223,107 +222,3 @@ def diff(self, variable):
             raise pybamm.ModelError("cannot take second time derivative of a Variable")
         else:
             return pybamm.Scalar(0)
-
-
-class ExternalVariable(Variable):
-    """
-    A node in the expression tree representing an external variable variable.
-
-    This node will be discretised by :class:`.Discretisation` and converted
-    to a :class:`.Vector` node.
-
-    Parameters
-    ----------
-
-    name : str
-        name of the node
-    domain : iterable of str
-        list of domains that this variable is valid over
-    auxiliary_domains : dict
-        dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...,
-        'quaternary': ...}). For example, for the single particle model, the particle
-        concentration would be a Variable with domain 'negative particle' and secondary
-        auxiliary domain 'current collector'. For the DFN, the particle concentration
-        would be a Variable with domain 'negative particle', secondary domain
-        'negative electrode' and tertiary domain 'current collector'
-    domains : dict
-        A dictionary equivalent to {'primary': domain, auxiliary_domains}. Either
-        'domain' and 'auxiliary_domains', or just 'domains', should be provided
-        (not both). In future, the 'domain' and 'auxiliary_domains' arguments may be
-        deprecated.
-    scale : float or :class:`pybamm.Symbol`, optional
-        The scale of the variable, used for scaling the model when solving. The state
-        vector representing this variable will be multiplied by this scale.
-        Default is 1.
-    reference : float or :class:`pybamm.Symbol`, optional
-        The reference value of the variable, used for scaling the model when solving.
-        This value will be added to the state vector representing this variable.
-        Default is 0.
-
-    *Extends:* :class:`pybamm.Variable`
-    """
-
-    def __init__(
-        self,
-        name,
-        size,
-        domain=None,
-        auxiliary_domains=None,
-        domains=None,
-        scale=1,
-        reference=0,
-    ):
-        self._size = size
-        super().__init__(
-            name, domain, auxiliary_domains, domains, scale=scale, reference=reference
-        )
-
-    @property
-    def size(self):
-        return self._size
-
-    def create_copy(self):
-        """See :meth:`pybamm.Symbol.new_copy()`."""
-        return ExternalVariable(self.name, self.size, domains=self.domains)
-
-    def _evaluate_for_shape(self):
-        """See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()`"""
-        return np.nan * np.ones((self.size, 1))
-
-    def _base_evaluate(self, t=None, y=None, y_dot=None, inputs=None):
-        # inputs should be a dictionary
-        # convert 'None' to empty dictionary for more informative error
-        if inputs is None:
-            inputs = {}
-        if not isinstance(inputs, dict):
-            # if the special input "shape test" is passed, just return 1
-            if inputs == "shape test":
-                return self.evaluate_for_shape()
-            raise TypeError("inputs should be a dictionary")
-        try:
-            out = inputs[self.name]
-            if isinstance(out, numbers.Number) or out.shape[0] == 1:
-                return out * np.ones((self.size, 1))
-            elif out.shape[0] != self.size:
-                raise ValueError(
-                    "External variable input has size {} but should be {}".format(
-                        out.shape[0], self.size
-                    )
-                )
-            else:
-                if isinstance(out, np.ndarray) and out.ndim == 1:
-                    out = out[:, np.newaxis]
-                return out
-        # raise more informative error if can't find name in dict
-        except KeyError:
-            raise KeyError("External variable '{}' not found".format(self.name))
-
-    def diff(self, variable):
-        if variable == self:
-            return pybamm.Scalar(1)
-        elif variable == pybamm.t:
-            raise pybamm.ModelError(
-                "cannot take time derivative of an external variable"
-            )
-        else:
-            return pybamm.Scalar(0)
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 31abb50cda..38b957f40f 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -90,9 +90,9 @@ class BaseModel:
 
     def __init__(self, name="Unnamed model"):
         self.name = name
-        self._options = {"external submodels": []}
+        self._options = {}
         self._built = False
-        self._built_fundamental_and_external = False
+        self._built_fundamental = False
 
         # Initialise empty model
         self.submodels = {}
@@ -109,7 +109,6 @@ def __init__(self, name="Unnamed model"):
         self._mass_matrix_inv = None
         self._jacobian = None
         self._jacobian_algebraic = None
-        self.external_variables = []
         self._parameters = None
         self._input_parameters = None
         self._parameter_info = None
@@ -419,7 +418,6 @@ def new_copy(self):
         new_model._boundary_conditions = self.boundary_conditions.copy()
         new_model._variables = self.variables.copy()
         new_model._events = self.events.copy()
-        new_model.external_variables = self.external_variables.copy()
         new_model._variables_casadi = self._variables_casadi.copy()
         return new_model
 
@@ -445,7 +443,7 @@ def update(self, *submodels):
             self.variables.update(submodel.variables)  # keys are strings so no check
             self._events += submodel.events
 
-    def build_fundamental_and_external(self):
+    def build_fundamental(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -455,23 +453,7 @@ def build_fundamental_and_external(self):
             )
             self.variables.update(submodel.get_fundamental_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
+        self._built_fundamental = True
 
     def build_coupled_variables(self):
         # Note: pybamm will try to get the coupled variables for the submodels in the
@@ -520,38 +502,37 @@ def build_coupled_variables(self):
     def build_model_equations(self):
         # Set model equations
         for submodel_name, submodel in self.submodels.items():
-            if submodel.external is False:
-                pybamm.logger.verbose(
-                    "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
-                )
+            pybamm.logger.verbose(
+                "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
+            )
 
-                submodel.set_rhs(self.variables)
-                pybamm.logger.verbose(
-                    "Setting algebraic for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_rhs(self.variables)
+            pybamm.logger.verbose(
+                "Setting algebraic for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_algebraic(self.variables)
-                pybamm.logger.verbose(
-                    "Setting boundary conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_algebraic(self.variables)
+            pybamm.logger.verbose(
+                "Setting boundary conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_boundary_conditions(self.variables)
-                pybamm.logger.verbose(
-                    "Setting initial conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
-                )
-                submodel.set_initial_conditions(self.variables)
-                submodel.set_events(self.variables)
-                pybamm.logger.verbose(
-                    "Updating {} submodel ({})".format(submodel_name, self.name)
+            submodel.set_boundary_conditions(self.variables)
+            pybamm.logger.verbose(
+                "Setting initial conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
-                self.update(submodel)
-                self.check_no_repeated_keys()
+            )
+            submodel.set_initial_conditions(self.variables)
+            submodel.set_events(self.variables)
+            pybamm.logger.verbose(
+                "Updating {} submodel ({})".format(submodel_name, self.name)
+            )
+            self.update(submodel)
+            self.check_no_repeated_keys()
 
     def build_model(self):
         self._build_model()
@@ -568,8 +549,8 @@ def _build_model(self):
 
         pybamm.logger.info("Start building {}".format(self.name))
 
-        if self._built_fundamental_and_external is False:
-            self.build_fundamental_and_external()
+        if self._built_fundamental is False:
+            self.build_fundamental()
 
         self.build_coupled_variables()
 
@@ -813,16 +794,7 @@ def check_well_determined(self, post_discretisation):
         all_vars_in_keys = all_vars_in_rhs_keys.union(all_vars_in_algebraic_keys)
         extra_variables_in_equations = all_vars_in_eqns.difference(all_vars_in_keys)
 
-        # get external variables
-        external_vars = set(self.external_variables)
-        for var in self.external_variables:
-            if isinstance(var, pybamm.Concatenation):
-                child_vars = set(var.children)
-                external_vars = external_vars.union(child_vars)
-
-        extra_variables = extra_variables_in_equations.difference(external_vars)
-
-        if extra_variables:
+        if extra_variables_in_equations:
             raise pybamm.ModelError("model is underdetermined (too many variables)")
 
     def check_algebraic_equations(self, post_discretisation):
@@ -858,13 +830,9 @@ def check_variables(self):
 
         vars_in_keys = set()
 
-        model_and_external_variables = (
-            list(self.rhs.keys())
-            + list(self.algebraic.keys())
-            + self.external_variables
-        )
+        model_keys = list(self.rhs.keys()) + list(self.algebraic.keys())
 
-        for var in model_and_external_variables:
+        for var in model_keys:
             if isinstance(var, pybamm.Variable):
                 vars_in_keys.add(var)
             # Key can be a concatenation
@@ -876,8 +844,8 @@ def check_variables(self):
                 raise pybamm.ModelError(
                     """
                     No key set for variable '{}'. Make sure it is included in either
-                    model.rhs, model.algebraic, or model.external_variables in an
-                    unmodified form (e.g. not Broadcasted)
+                    model.rhs or model.algebraic, in an unmodified form
+                    (e.g. not Broadcasted)
                     """.format(
                         var
                     )
@@ -948,8 +916,10 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         variable_names : list
             Variables to be exported alongside the model structure
         input_parameter_order : list, optional
-            Order in which the input parameters should be stacked. If None, the order
-            returned by :meth:`BaseModel.input_parameters` is used
+            Order in which the input parameters should be stacked.
+            If input_parameter_order=None and len(self.input_parameters) > 1, a
+            ValueError is raised (this helps to avoid accidentally using the wrong
+            order)
 
         Returns
         -------
@@ -970,21 +940,18 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         for input_param in self.input_parameters:
             name = input_param.name
             inputs_wrong_order[name] = casadi.MX.sym(name, input_param._expected_size)
-        # Read external variables
-        external_casadi = {}
-        for external_varaiable in self.external_variables:
-            name = external_varaiable.name
-            ev_size = external_varaiable._evaluate_for_shape().shape[0]
-            external_casadi[name] = casadi.MX.sym(name, ev_size)
         # Sort according to input_parameter_order
         if input_parameter_order is None:
+            if len(inputs_wrong_order) > 1:
+                raise ValueError(
+                    "input_parameter_order must be specified if there is more than one "
+                    "input parameter"
+                )
             inputs = inputs_wrong_order
         else:
             inputs = {name: inputs_wrong_order[name] for name in input_parameter_order}
-        # Set up external variables and inputs
-        # Put external variables first like the integrator expects
-        ext_and_in = {**external_casadi, **inputs}
-        inputs_stacked = casadi.vertcat(*[p for p in ext_and_in.values()])
+        # Set up inputs
+        inputs_stacked = casadi.vertcat(*[p for p in inputs.values()])
 
         # Convert initial conditions to casadi form
         y0 = self.concatenated_initial_conditions.to_casadi(
@@ -994,7 +961,7 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         z0 = y0[self.concatenated_rhs.size :]
 
         # Convert rhs and algebraic to casadi form and calculate jacobians
-        rhs = self.concatenated_rhs.to_casadi(t_casadi, y_casadi, inputs=ext_and_in)
+        rhs = self.concatenated_rhs.to_casadi(t_casadi, y_casadi, inputs=inputs)
         jac_rhs = casadi.jacobian(rhs, y_casadi)
         algebraic = self.concatenated_algebraic.to_casadi(
             t_casadi, y_casadi, inputs=inputs
@@ -1005,7 +972,7 @@ def export_casadi_objects(self, variable_names, input_parameter_order=None):
         variables = OrderedDict()
         for name in variable_names:
             var = self.variables[name]
-            variables[name] = var.to_casadi(t_casadi, y_casadi, inputs=ext_and_in)
+            variables[name] = var.to_casadi(t_casadi, y_casadi, inputs=inputs)
 
         casadi_dict = {
             "t": t_casadi,
@@ -1036,8 +1003,10 @@ def generate(
         variable_names : list
             Variables to be exported alongside the model structure
         input_parameter_order : list, optional
-            Order in which the input parameters should be stacked. If None, the order
-            returned by :meth:`BaseModel.input_parameters` is used
+            Order in which the input parameters should be stacked.
+            If input_parameter_order=None and len(self.input_parameters) > 1, a
+            ValueError is raised (this helps to avoid accidentally using the wrong
+            order)
         cg_options : dict
             Options to pass to the code generator.
             See https://web.casadi.org/docs/#generating-c-code
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index ac254ce892..b9d177c8ee 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -50,11 +50,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "electrolyte conductivity" : str
                 Can be "default" (default), "full", "leading order", "composite" or
                 "integrated".
-            * "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`.
             * "hydrolysis" : str
                 Whether to include hydrolysis in the model. Only implemented for
                 lead-acid models. Can be "false" (default) or "true". If "true", then
@@ -278,7 +273,6 @@ def __init__(self, extra_options):
         default_options = {
             name: options[0] for name, options in self.possible_options.items()
         }
-        default_options["external submodels"] = []
         default_options["timescale"] = "default"
 
         # Change the default for cell geometry based on which thermal option is provided
@@ -530,7 +524,7 @@ def __init__(self, extra_options):
 
         # Check options are valid
         for option, value in options.items():
-            if option in ["external submodels", "working electrode"]:
+            if option in ["working electrode"]:
                 pass
             else:
                 if isinstance(value, str) or option in [
@@ -890,7 +884,7 @@ def set_standard_output_variables(self):
                 {"y": var.y, "y [m]": var.y * L_z, "z": var.z, "z [m]": var.z * L_z}
             )
 
-    def build_fundamental_and_external(self):
+    def build_fundamental(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
             pybamm.logger.debug(
@@ -900,23 +894,7 @@ def build_fundamental_and_external(self):
             )
             self.variables.update(submodel.get_fundamental_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
+        self._built_fundamental = True
 
     def build_coupled_variables(self):
         # Note: pybamm will try to get the coupled variables for the submodels in the
@@ -965,38 +943,37 @@ def build_coupled_variables(self):
     def build_model_equations(self):
         # Set model equations
         for submodel_name, submodel in self.submodels.items():
-            if submodel.external is False:
-                pybamm.logger.verbose(
-                    "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
-                )
+            pybamm.logger.verbose(
+                "Setting rhs for {} submodel ({})".format(submodel_name, self.name)
+            )
 
-                submodel.set_rhs(self.variables)
-                pybamm.logger.verbose(
-                    "Setting algebraic for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_rhs(self.variables)
+            pybamm.logger.verbose(
+                "Setting algebraic for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_algebraic(self.variables)
-                pybamm.logger.verbose(
-                    "Setting boundary conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
+            submodel.set_algebraic(self.variables)
+            pybamm.logger.verbose(
+                "Setting boundary conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
+            )
 
-                submodel.set_boundary_conditions(self.variables)
-                pybamm.logger.verbose(
-                    "Setting initial conditions for {} submodel ({})".format(
-                        submodel_name, self.name
-                    )
-                )
-                submodel.set_initial_conditions(self.variables)
-                submodel.set_events(self.variables)
-                pybamm.logger.verbose(
-                    "Updating {} submodel ({})".format(submodel_name, self.name)
+            submodel.set_boundary_conditions(self.variables)
+            pybamm.logger.verbose(
+                "Setting initial conditions for {} submodel ({})".format(
+                    submodel_name, self.name
                 )
-                self.update(submodel)
-                self.check_no_repeated_keys()
+            )
+            submodel.set_initial_conditions(self.variables)
+            submodel.set_events(self.variables)
+            pybamm.logger.verbose(
+                "Updating {} submodel ({})".format(submodel_name, self.name)
+            )
+            self.update(submodel)
+            self.check_no_repeated_keys()
 
     def build_model(self):
 
diff --git a/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py b/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
index 078bb5de29..0414f3fde3 100644
--- a/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
+++ b/pybamm/models/full_battery_models/equivalent_circuit/thevenin.py
@@ -44,11 +44,6 @@ class Thevenin(pybamm.BaseModel):
                 - callable : if a callable is given as this option, the function \
                     defines the residual of an algebraic equation. The applied current \
                     will be solved for such that the algebraic constraint is satisfied.
-            * "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`.
     build :  bool, optional
         Whether to build the model on instantiation. Default is True. Setting this
         option to False allows users to change any number of the submodels before
@@ -89,7 +84,6 @@ def set_options(self, extra_options=None):
             "calculate discharge energy": ["false", "true"],
             "operating mode": OperatingModes("current"),
             "number of rc elements": NaturalNumberOption(1),
-            "external submodels": [[]],
         }
 
         default_options = {
diff --git a/pybamm/models/submodels/active_material/constant_active_material.py b/pybamm/models/submodels/active_material/constant_active_material.py
index 66219ee4ab..2943cbbfd1 100644
--- a/pybamm/models/submodels/active_material/constant_active_material.py
+++ b/pybamm/models/submodels/active_material/constant_active_material.py
@@ -37,7 +37,7 @@ def get_fundamental_variables(self):
 
         variables.update(
             {
-                "Loss of lithium due to LAM "
+                "Loss of lithium due to loss of active material "
                 f"in {domain} electrode [mol]": pybamm.Scalar(0)
             }
         )
diff --git a/pybamm/models/submodels/active_material/loss_active_material.py b/pybamm/models/submodels/active_material/loss_active_material.py
index 9c68aeb96d..fce67148f6 100644
--- a/pybamm/models/submodels/active_material/loss_active_material.py
+++ b/pybamm/models/submodels/active_material/loss_active_material.py
@@ -54,11 +54,12 @@ def get_fundamental_variables(self):
             )
         variables = self._get_standard_active_material_variables(eps_solid)
         lli_due_to_lam = pybamm.Variable(
-            f"Loss of lithium due to LAM in {domain} electrode [mol]"
+            "Loss of lithium due to loss of active material "
+            f"in {domain} electrode [mol]"
         )
         variables.update(
             {
-                "Loss of lithium due to LAM "
+                "Loss of lithium due to loss of active material "
                 f"in {domain} electrode [mol]": lli_due_to_lam
             }
         )
@@ -142,12 +143,13 @@ def set_rhs(self, variables):
                 f"{Domain} electrode active material volume fraction change"
             ]
 
-        # Loss of lithium due to LAM
+        # Loss of lithium due to loss of active material
         # See eq 37 in "Sulzer, Valentin, et al. "Accelerated battery lifetime
         # simulations using adaptive inter-cycle extrapolation algorithm."
         # Journal of The Electrochemical Society 168.12 (2021): 120531.
         lli_due_to_lam = variables[
-            f"Loss of lithium due to LAM in {domain} electrode [mol]"
+            "Loss of lithium due to loss of active material "
+            f"in {domain} electrode [mol]"
         ]
         # Multiply by mol.m-3 * m3 to get mol
         c_s_av = variables[f"Average {domain} particle concentration [mol.m-3]"]
@@ -174,6 +176,7 @@ def set_initial_conditions(self, variables):
             self.initial_conditions = {eps_solid: eps_solid_init}
 
         lli_due_to_lam = variables[
-            f"Loss of lithium due to LAM in {domain} electrode [mol]"
+            "Loss of lithium due to loss of active material "
+            f"in {domain} electrode [mol]"
         ]
         self.initial_conditions[lli_due_to_lam] = pybamm.Scalar(0)
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index bbabc0d4da..a2f450e04b 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -150,52 +150,6 @@ def get_fundamental_variables(self):
         """
         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]
-                    # also add the ids of the children to unique ids
-                    for child in var.children:
-                        unique_ids += [child]
-
-            # 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]
-
-        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/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index faba47b98d..9454f1ab82 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -470,11 +470,6 @@ def process_model(self, unprocessed_model, inplace=True):
 
         model.events = new_events
 
-        # Set external variables
-        model.external_variables = [
-            self.process_symbol(var) for var in unprocessed_model.external_variables
-        ]
-
         # Process timescale
         new_timescale = self.process_symbol(unprocessed_model.timescale)
         if isinstance(new_timescale, pybamm.Scalar):
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 972debf03f..b41b48581b 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -108,7 +108,7 @@ def __init__(
         self.submesh_types = submesh_types or self.model.default_submesh_types
         self.var_pts = var_pts or self.model.default_var_pts
         self.spatial_methods = spatial_methods or self.model.default_spatial_methods
-        self._solver = solver or self.model.default_solver
+        self.solver = solver or self.model.default_solver
         self.output_variables = output_variables
 
         # Initialize empty built states
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 6e85d067ab..7c4b20cd75 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -654,7 +654,6 @@ def solve(
         self,
         model,
         t_eval=None,
-        external_variables=None,
         inputs=None,
         initial_conditions=None,
         nproc=None,
@@ -672,9 +671,6 @@ def solve(
             an error is raised
         t_eval : numeric type
             The times (in seconds) at which to compute the solution
-        external_variables : dict
-            A dictionary of external variables and their corresponding
-            values at the current time
         inputs : dict or list, optional
             A dictionary or list of dictionaries describing any input parameters to
             pass to the model when solving
@@ -749,16 +745,15 @@ def solve(
         if (np.diff(t_eval) < 0).any():
             raise pybamm.SolverError("t_eval must increase monotonically")
 
-        # Set up external variables and inputs
+        # Set up inputs
         #
         # Argument "inputs" can be either a list of input dicts or
         # a single dict. The remaining of this function is only working
         # with variable "input_list", which is a list of dictionaries.
         # If "inputs" is a single dict, "inputs_list" is a list of only one dict.
         inputs_list = inputs if isinstance(inputs, list) else [inputs]
-        ext_and_inputs_list = [
-            self._set_up_ext_and_inputs(model, external_variables, inputs)
-            for inputs in inputs_list
+        model_inputs_list = [
+            self._set_up_model_inputs(model, inputs) for inputs in inputs_list
         ]
 
         # Cannot use multiprocessing with model in "jax" format
@@ -795,10 +790,10 @@ def solve(
                 )
             # It is assumed that when len(inputs_list) > 1, model set
             # up (initial condition, time-scale and length-scale) does
-            # not depend on input parameters. Thefore only `ext_and_inputs[0]`
+            # not depend on input parameters. Thefore only `model_inputs[0]`
             # is passed to `set_up`.
             # See https://github.com/pybamm-team/PyBaMM/pull/1261
-            self.set_up(model, ext_and_inputs_list[0], t_eval)
+            self.set_up(model, model_inputs_list[0], t_eval)
             self._model_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
@@ -814,7 +809,7 @@ def solve(
                 else:
                     # If the new initial conditions are different
                     # and cannot be evaluated directly, set up again
-                    self.set_up(model, ext_and_inputs_list[0], t_eval, ics_only=True)
+                    self.set_up(model, model_inputs_list[0], t_eval, ics_only=True)
                 self._model_set_up[model][
                     "initial conditions"
                 ] = model.concatenated_initial_conditions
@@ -824,13 +819,13 @@ def solve(
 
         # (Re-)calculate consistent initial conditions
         # Assuming initial conditions do not depend on input parameters
-        # when len(inputs_list) > 1, only `ext_and_inputs_list[0]`
+        # when len(inputs_list) > 1, only `model_inputs_list[0]`
         # is passed to `_set_initial_conditions`.
         # See https://github.com/pybamm-team/PyBaMM/pull/1261
         if len(inputs_list) > 1:
             all_inputs_names = set(
                 itertools.chain.from_iterable(
-                    [ext_and_inputs.keys() for ext_and_inputs in ext_and_inputs_list]
+                    [model_inputs.keys() for model_inputs in model_inputs_list]
                 )
             )
             initial_conditions_node_names = set(
@@ -846,12 +841,12 @@ def solve(
         t_eval_dimensionless = t_eval / model.timescale_eval
 
         self._set_initial_conditions(
-            model, t_eval_dimensionless[0], ext_and_inputs_list[0], update_rhs=True
+            model, t_eval_dimensionless[0], model_inputs_list[0], update_rhs=True
         )
 
         # Check initial conditions don't violate events
         self._check_events_with_initial_conditions(
-            t_eval_dimensionless, model, ext_and_inputs_list[0]
+            t_eval_dimensionless, model, model_inputs_list[0]
         )
 
         # Process discontinuities
@@ -875,12 +870,12 @@ def solve(
                     t_eval_dimensionless[end_index - 1] * model.timescale_eval,
                 )
             )
-            ninputs = len(ext_and_inputs_list)
+            ninputs = len(model_inputs_list)
             if ninputs == 1:
                 new_solution = self._integrate(
                     model,
                     t_eval_dimensionless[start_index:end_index],
-                    ext_and_inputs_list[0],
+                    model_inputs_list[0],
                 )
                 new_solutions = [new_solution]
             else:
@@ -890,7 +885,7 @@ def solve(
                         zip(
                             [model] * ninputs,
                             [t_eval_dimensionless[start_index:end_index]] * ninputs,
-                            ext_and_inputs_list,
+                            model_inputs_list,
                         ),
                     )
                     p.close()
@@ -917,7 +912,7 @@ def solve(
                 model.y0 = last_state
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], ext_and_inputs_list[0]
+                        model, t_eval_dimensionless[end_index], model_inputs_list[0]
                     )
         solve_time = timer.time()
 
@@ -1077,7 +1072,6 @@ def step(
         model,
         dt,
         npts=2,
-        external_variables=None,
         inputs=None,
         save=True,
     ):
@@ -1098,9 +1092,6 @@ def step(
         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
         inputs : dict, optional
             Any input parameters to pass to the model when solving
         save : bool
@@ -1142,8 +1133,8 @@ def step(
         # Set timer
         timer = pybamm.Timer()
 
-        # Set up external variables and inputs
-        ext_and_inputs = self._set_up_ext_and_inputs(model, external_variables, inputs)
+        # Set up inputs
+        model_inputs = self._set_up_model_inputs(model, inputs)
 
         t = old_solution.t[-1]
 
@@ -1157,7 +1148,7 @@ def step(
                     f'"{existing_model.name}". Please create a separate '
                     "solver for this model"
                 )
-            self.set_up(model, ext_and_inputs)
+            self.set_up(model, model_inputs)
             self._model_set_up.update(
                 {model: {"initial conditions": model.concatenated_initial_conditions}}
             )
@@ -1173,7 +1164,7 @@ def step(
         if isinstance(old_solution, pybamm.EmptySolution):
             if not first_step_this_model:
                 # reset y0 to original initial conditions
-                self.set_up(model, ext_and_inputs, ics_only=True)
+                self.set_up(model, model_inputs, ics_only=True)
         else:
             if old_solution.all_models[-1] == model:
                 # initialize with old solution
@@ -1183,20 +1174,20 @@ def step(
                     old_solution, return_type="ics"
                 )
                 model.y0 = concatenated_initial_conditions.evaluate(
-                    0, inputs=ext_and_inputs
+                    0, inputs=model_inputs
                 )
 
         set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions
-        self._set_initial_conditions(model, t, ext_and_inputs, update_rhs=False)
+        self._set_initial_conditions(model, t, model_inputs, update_rhs=False)
 
         # Non-dimensionalise dt
         dt_dimensionless = dt / model.timescale_eval
         t_eval = np.linspace(t, t + dt_dimensionless, npts)
 
         # Check initial conditions don't violate events
-        self._check_events_with_initial_conditions(t_eval, model, ext_and_inputs)
+        self._check_events_with_initial_conditions(t_eval, model, model_inputs)
 
         # Step
         pybamm.logger.verbose(
@@ -1206,7 +1197,7 @@ def step(
             )
         )
         timer.reset()
-        solution = self._integrate(model, t_eval, ext_and_inputs)
+        solution = self._integrate(model, t_eval, model_inputs)
         solution.solve_time = timer.time()
 
         # Check if extrapolation occurred
@@ -1364,8 +1355,8 @@ def check_extrapolation(self, solution, events):
                         "outside these bounds."
                     )
 
-    def _set_up_ext_and_inputs(self, model, external_variables, inputs):
-        """Set up external variables and input parameters"""
+    def _set_up_model_inputs(self, model, inputs):
+        """Set up input parameters"""
         inputs = inputs or {}
 
         # Go through all input parameters that can be found in the model
@@ -1380,14 +1371,11 @@ def _set_up_ext_and_inputs(self, model, external_variables, inputs):
                 raise pybamm.SolverError(f"No value provided for input '{name}'")
         inputs = inputs_in_model
 
-        external_variables = external_variables or {}
-
         ordered_inputs_names = list(inputs.keys())
         ordered_inputs_names.sort()
         ordered_inputs = {name: inputs[name] for name in ordered_inputs_names}
 
-        ext_and_inputs = {**external_variables, **ordered_inputs}
-        return ext_and_inputs
+        return ordered_inputs
 
 
 def process(symbol, name, vars_for_processing, use_jacobian=None):
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index eea853eeda..12627b7120 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -143,7 +143,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         t_eval : numeric type
             The times at which to compute the solution
         inputs_dict : dict, optional
-            Any external variables or input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving
         """
 
         # Record whether there are any symbolic inputs
@@ -605,7 +605,7 @@ def _run_integrator(
         y0:
             casadi vector of initial conditions
         inputs_dict : dict, optional
-            Any external variables or input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving
         inputs:
             Casadi vector of inputs
         t_eval : numeric type
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 8a7671f084..1c2822d968 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -510,7 +510,7 @@ def _integrate(self, model, t_eval, inputs_dict=None):
         t_eval : numeric type
             The times at which to compute the solution
         inputs_dict : dict, optional
-            Any external variables or input parameters to pass to the model when solving
+            Any input parameters to pass to the model when solving
         """
         inputs_dict = inputs_dict or {}
         # stack inputs
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 8b069ff626..423cdbe356 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -65,32 +65,6 @@ def spatial_variable(self, symbol):
             entries = np.tile(symbol_mesh.nodes, repeats)
         return pybamm.Vector(entries, domains=symbol.domains)
 
-    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(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the gradient operator.
         See :meth:`pybamm.SpatialMethod.gradient`
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 6826414bab..844fb9559e 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -330,9 +330,6 @@ def internal_neumann_condition(
 
         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/pybamm/spatial_methods/spectral_volume.py b/pybamm/spatial_methods/spectral_volume.py
index 17fe70e040..af62570cc8 100644
--- a/pybamm/spatial_methods/spectral_volume.py
+++ b/pybamm/spatial_methods/spectral_volume.py
@@ -14,7 +14,7 @@ class SpectralVolume(pybamm.FiniteVolume):
     subdivision of any 1D mesh, so it shouldn't be a problem).
 
     For broadcast and mass_matrix, we follow the default behaviour from
-    SpatialMethod. For spatial_variable, preprocess_external_variables,
+    SpatialMethod. For spatial_variable,
     divergence, divergence_matrix, laplacian, integral,
     definite_integral_matrix, indefinite_integral,
     indefinite_integral_matrix, indefinite_integral_matrix_nodes,
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 4621cfecc5..8dd6cccbe7 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -295,10 +295,10 @@ def __init__(self, model, param, disc, solution, operating_condition):
 
         self.n_Li_side = solution["Total lithium lost to side reactions [mol]"]
         self.n_Li_LAM_n = solution[
-            "Loss of lithium due to LAM in negative electrode [mol]"
+            "Loss of lithium due to loss of active material in negative electrode [mol]"
         ]
         self.n_Li_LAM_p = solution[
-            "Loss of lithium due to LAM in positive electrode [mol]"
+            "Loss of lithium due to loss of active material in positive electrode [mol]"
         ]
 
         if model.options["particle size"] == "distribution":
@@ -807,10 +807,10 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.n_Li_lost = solution["Total lithium lost [mol]"]
         self.n_Li_lost_rxn = solution["Total lithium lost to side reactions [mol]"]
         self.n_Li_lost_LAM_n = solution[
-            "Loss of lithium due to LAM in negative electrode [mol]"
+            "Loss of lithium due to loss of active material in negative electrode [mol]"
         ]
         self.n_Li_lost_LAM_p = solution[
-            "Loss of lithium due to LAM in positive electrode [mol]"
+            "Loss of lithium due to loss of active material in positive electrode [mol]"
         ]
 
     def test_degradation_modes(self):
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
deleted file mode 100644
index 977c46a618..0000000000
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py
+++ /dev/null
@@ -1,57 +0,0 @@
-#
-# Tests for external submodels
-#
-import pybamm
-import unittest
-import numpy as np
-
-
-class TestExternalCC(unittest.TestCase):
-    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,
-        }
-        sim = pybamm.Simulation(model, var_pts=var_pts)
-
-        # Simulate 100 seconds
-        t_eval = np.linspace(0, 100, 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.solution["Positive current collector potential"].data[:, -1]
-
-        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
index f6b00d0564..9a1c8ff732 100644
--- 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
@@ -1,92 +1,31 @@
 #
-# Tests for external submodels
+# Tests for inputting a temperature profile
 #
 import pybamm
 import unittest
 import numpy as np
 
 
-class TestExternalThermalModels(unittest.TestCase):
-    def test_external_lumped_temperature(self):
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
+class TestInputLumpedTemperature(unittest.TestCase):
+    def test_input_lumped_temperature(self):
+        model = pybamm.lithium_ion.SPMe()
+        parameter_values = model.default_parameter_values
+        # in the default isothermal model, the temperature is everywhere equal
+        # to the ambient temperature
+        parameter_values["Ambient temperature [K]"] = pybamm.InputParameter(
+            "Volume-averaged cell temperature [K]"
+        )
         sim = pybamm.Simulation(model)
 
         t_eval = np.linspace(0, 100, 3)
 
-        T_av = 0
+        T_av = 298
 
         for i in np.arange(1, len(t_eval) - 1):
             dt = t_eval[i + 1] - t_eval[i]
-            external_variables = {"Volume-averaged cell temperature": T_av}
+            inputs = {"Volume-averaged cell temperature [K]": T_av}
             T_av += 1
-            sim.step(dt, external_variables=external_variables)
-
-    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, 100, 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)
-
-    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.CasadiSolver()
-        sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver)
-        sim.build()
-
-        t_eval = np.linspace(0, 100, 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)
+            sim.step(dt, inputs=inputs)  # works
 
 
 if __name__ == "__main__":
diff --git a/tests/integration/test_solvers/test_external_variables.py b/tests/integration/test_solvers/test_external_variables.py
deleted file mode 100644
index 130023f3b4..0000000000
--- a/tests/integration/test_solvers/test_external_variables.py
+++ /dev/null
@@ -1,59 +0,0 @@
-#
-# Test solvers with external variables
-#
-import pybamm
-import numpy as np
-import os
-import unittest
-
-
-class TestExternalVariables(unittest.TestCase):
-    def test_on_dfn(self):
-        e_height = 0.25
-
-        model = pybamm.lithium_ion.DFN()
-        geometry = model.default_geometry
-        param = model.default_parameter_values
-        param.update({"Negative electrode conductivity [S.m-1]": "[input]"})
-        inputs = {"Negative electrode conductivity [S.m-1]": e_height}
-        var_pts = {"x_n": 5, "x_s": 5, "x_p": 5, "r_n": 10, "r_p": 10}
-        spatial_methods = model.default_spatial_methods
-
-        solver = pybamm.CasadiSolver()
-        sim = pybamm.Simulation(
-            model=model,
-            geometry=geometry,
-            parameter_values=param,
-            var_pts=var_pts,
-            spatial_methods=spatial_methods,
-            solver=solver,
-        )
-        sim.solve(t_eval=np.linspace(0, 3600, 100), inputs=inputs)
-
-    def test_external_variables_SPMe_thermal(self):
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-        t_eval = np.linspace(0, 100, 3)
-        T_av = 0
-        for i in np.arange(1, len(t_eval) - 1):
-            dt = t_eval[i + 1] - t_eval[i]
-            external_variables = {"Volume-averaged cell temperature": T_av}
-            T_av += 1
-            sim.step(dt, external_variables=external_variables)
-        V = sim.solution["Terminal voltage [V]"].data
-        np.testing.assert_array_less(np.diff(V), 0)
-        # test generate with external variable
-        sim.built_model.generate("test.c", ["Volume-averaged cell temperature"])
-        os.remove("test.c")
-
-
-if __name__ == "__main__":
-    import sys
-
-    print("Add -v for more debug output")
-
-    if "-v" in sys.argv:
-        debug = True
-    pybamm.settings.debug_mode = True
-    unittest.main()
diff --git a/tests/integration/test_solvers/test_solution.py b/tests/integration/test_solvers/test_solution.py
index 3493f66e46..f2837e1a7b 100644
--- a/tests/integration/test_solvers/test_solution.py
+++ b/tests/integration/test_solvers/test_solution.py
@@ -50,64 +50,6 @@ def test_append(self):
             decimal=4,
         )
 
-    def test_append_external_variables(self):
-        model = pybamm.lithium_ion.SPM(
-            {
-                "thermal": "lumped",
-                "external submodels": ["thermal", "negative primary particle"],
-            }
-        )
-        # 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
-        mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
-
-        # discretise model
-        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-        disc.process_model(model)
-
-        # solve model
-        solver = model.default_solver
-
-        Nr = model.default_var_pts["r_n"]
-
-        T_av = 0
-        c_s_n_av = np.ones((Nr, 1)) * 0.6
-        external_variables = {
-            "Volume-averaged cell temperature": T_av,
-            "X-averaged negative particle concentration": c_s_n_av,
-        }
-
-        # Step
-        dt = 0.1
-        sol_step = None
-        for _ in range(5):
-            sol_step = solver.step(
-                sol_step, model, dt, external_variables=external_variables
-            )
-        np.testing.assert_array_equal(
-            sol_step.all_inputs[0]["Volume-averaged cell temperature"], 0
-        )
-        np.testing.assert_array_equal(
-            sol_step.all_inputs[0]["X-averaged negative particle concentration"], 0.6
-        )
-
-        # Solve
-        t_eval = np.linspace(0, 3600)
-        sol = solver.solve(model, t_eval, external_variables=external_variables)
-        np.testing.assert_array_equal(
-            sol.all_inputs[0]["Volume-averaged cell temperature"], 0
-        )
-        np.testing.assert_array_equal(
-            sol.all_inputs[0]["X-averaged negative particle concentration"], 0.6
-        )
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 01dd02a2b4..295d8a2780 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -67,184 +67,6 @@ def test_add_internal_boundary_conditions(self):
         for child in c_e.children:
             self.assertTrue(child 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.assertIsInstance(model.variables["b"], pybamm.ExternalVariable)
-        self.assertEqual(model.variables["b"].evaluate(inputs={"b": np.array([1])}), 1)
-
-    def test_adding_0D_external_variable_fail(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]
-
-        disc = pybamm.Discretisation()
-        with self.assertRaisesRegex(ValueError, "Variable b must be in the model"):
-            disc.process_model(model)
-
-    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": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
-
-        submesh_types = {"test": 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(disc.y_slices[a][0], slice(0, 10, None))
-
-        self.assertEqual(model.y_slices[a][0], slice(0, 10, None))
-        self.assertEqual(model.bounds, disc.bounds)
-
-        b_test = np.ones((10, 1))
-        np.testing.assert_array_equal(
-            model.variables["b"].evaluate(inputs={"b": b_test}), b_test
-        )
-
-        # check that b is added to the boundary conditions
-        model.bcs[b]["left"]
-        model.bcs[b]["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_concatenation_external_variables(self):
-        model = pybamm.BaseModel()
-
-        a = pybamm.Variable("a", domain=["test", "test1"])
-        b1 = pybamm.Variable("b1", domain=["test"])
-        b2 = pybamm.Variable("b2", domain=["test1"])
-        b = pybamm.concatenation(b1, b2)
-
-        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,
-            "b1": b1,
-            "b2": b2,
-            "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")
-        y = pybamm.SpatialVariable("y", domain="test1", coord_sys="cartesian")
-        geometry = {
-            "test": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}},
-            "test1": {y: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}},
-        }
-
-        submesh_types = {
-            "test": pybamm.Uniform1DSubMesh,
-            "test1": pybamm.Uniform1DSubMesh,
-        }
-        var_pts = {x: 10, y: 5}
-        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-
-        spatial_methods = {
-            "test": pybamm.FiniteVolume(),
-            "test1": pybamm.FiniteVolume(),
-        }
-        disc = pybamm.Discretisation(mesh, spatial_methods)
-        disc.process_model(model)
-
-        self.assertEqual(disc.y_slices[a][0], slice(0, 15, None))
-
-        b_test = np.linspace(0, 1, 15)[:, np.newaxis]
-        np.testing.assert_array_equal(
-            model.variables["b"].evaluate(inputs={"b": b_test}), b_test
-        )
-        np.testing.assert_array_equal(
-            model.variables["b1"].evaluate(inputs={"b": b_test}), b_test[:10]
-        )
-        np.testing.assert_array_equal(
-            model.variables["b2"].evaluate(inputs={"b": b_test}), b_test[10:]
-        )
-
-        # check that b is added to the boundary conditions
-        model.bcs[b]["left"]
-        model.bcs[b]["right"]
-
-        # check that grad and div(grad ) produce the correct shapes
-        self.assertEqual(model.variables["b"].shape_for_testing, (15, 1))
-        self.assertEqual(model.variables["grad b"].shape_for_testing, (16, 1))
-        self.assertEqual(model.variables["div grad b"].shape_for_testing, (15, 1))
-        self.assertEqual(model.variables["b1"].shape_for_testing, (10, 1))
-        self.assertEqual(model.variables["b2"].shape_for_testing, (5, 1))
-
-    def test_adding_2D_external_variable_fail(self):
-        model = pybamm.BaseModel()
-        a = pybamm.Variable(
-            "a",
-            domain=["negative electrode", "separator"],
-            auxiliary_domains={"secondary": "current collector"},
-        )
-        b1 = pybamm.Variable(
-            "b",
-            domain="negative electrode",
-            auxiliary_domains={"secondary": "current collector"},
-        )
-        b2 = pybamm.Variable(
-            "b",
-            domain="separator",
-            auxiliary_domains={"secondary": "current collector"},
-        )
-        b = pybamm.concatenation(b1, b2)
-
-        model.rhs = {a: a * b}
-        model.initial_conditions = {a: 0}
-        model.external_variables = [b]
-        model.variables = {"b": b}
-
-        disc = get_1p1d_discretisation_for_testing()
-        with self.assertRaisesRegex(
-            NotImplementedError, "Cannot create 2D external variable"
-        ):
-            disc.process_model(model)
-
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()
diff --git a/tests/unit/test_expression_tree/test_d_dt.py b/tests/unit/test_expression_tree/test_d_dt.py
index ce4b1d920b..9171452406 100644
--- a/tests/unit/test_expression_tree/test_d_dt.py
+++ b/tests/unit/test_expression_tree/test_d_dt.py
@@ -40,9 +40,6 @@ def test_time_derivative_of_variable(self):
         with self.assertRaises(pybamm.ModelError):
             a = (pybamm.Variable("a")).diff(pybamm.t).diff(pybamm.t)
 
-        with self.assertRaises(pybamm.ModelError):
-            a = pybamm.ExternalVariable("a", 1).diff(pybamm.t)
-
     def test_time_derivative_of_state_vector(self):
 
         sv = pybamm.StateVector(slice(0, 10))
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 ba91ae8f08..7e3874eb07 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
@@ -360,31 +360,6 @@ def test_convert_input_parameter(self):
             casadi_inputs["Input 2"] * casadi_y,
         )
 
-    def test_convert_external_variable(self):
-        casadi_t = casadi.MX.sym("t")
-        casadi_y = casadi.MX.sym("y", 10)
-        casadi_inputs = {
-            "External 1": casadi.MX.sym("External 1", 3),
-            "External 2": casadi.MX.sym("External 2", 10),
-        }
-
-        pybamm_y = pybamm.StateVector(slice(0, 10))
-        pybamm_u1 = pybamm.ExternalVariable("External 1", 3)
-        pybamm_u2 = pybamm.ExternalVariable("External 2", 10)
-
-        # External only
-        self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, inputs=casadi_inputs),
-            casadi_inputs["External 1"],
-        )
-
-        # More complex
-        expr = pybamm_u2 + pybamm_y
-        self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, inputs=casadi_inputs),
-            casadi_inputs["External 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 684037b28a..d9558de3ae 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -55,12 +55,6 @@ def test_symbol_new_copy(self):
             pybamm.t,
             pybamm.Index(vec, 1),
             pybamm.NotConstant(a),
-            pybamm.ExternalVariable(
-                "external variable",
-                20,
-                domain="test",
-                auxiliary_domains={"secondary": "test2"},
-            ),
             pybamm.minimum(a, b),
             pybamm.maximum(a, b),
             pybamm.SparseStack(mat, mat),
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 1dbe3b67db..ee2b97e8d1 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -90,42 +90,6 @@ def test_variable_diff(self):
         self.assertEqual(a.diff(b).evaluate(), 0)
 
 
-class TestExternalVariable(unittest.TestCase):
-    def test_external_variable_scalar(self):
-        a = pybamm.ExternalVariable("a", 1)
-        self.assertEqual(a.size, 1)
-
-        self.assertEqual(a.evaluate(inputs={"a": 3}), 3)
-
-        with self.assertRaisesRegex(KeyError, "External variable"):
-            a.evaluate()
-        with self.assertRaisesRegex(TypeError, "inputs should be a dictionary"):
-            a.evaluate(inputs="not a dictionary")
-
-    def test_external_variable_vector(self):
-        a = pybamm.ExternalVariable("a", 10)
-        self.assertEqual(a.size, 10)
-
-        a_test = 2 * np.ones((10, 1))
-        np.testing.assert_array_equal(a.evaluate(inputs={"a": a_test}), a_test)
-        np.testing.assert_array_equal(
-            a.evaluate(inputs={"a": a_test.flatten()}), a_test
-        )
-
-        np.testing.assert_array_equal(a.evaluate(inputs={"a": 2}), a_test)
-
-        with self.assertRaisesRegex(ValueError, "External variable"):
-            a.evaluate(inputs={"a": np.ones((5, 1))})
-
-    def test_external_variable_diff(self):
-        a = pybamm.ExternalVariable("a", 10)
-        b = pybamm.Variable("b")
-        self.assertIsInstance(a.diff(a), pybamm.Scalar)
-        self.assertEqual(a.diff(a).evaluate(), 1)
-        self.assertIsInstance(a.diff(b), pybamm.Scalar)
-        self.assertEqual(a.diff(b).evaluate(), 0)
-
-
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index 161f347a26..1d00cca1b4 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -531,15 +531,33 @@ def test_export_casadi(self):
         np.testing.assert_array_equal(np.array(jac_alg_fn(5, 6, 7, [9, 8])), [[1, -1]])
         self.assertEqual(var_fn(6, 3, 2, [2, 7]), -1)
 
-        # Test model with external variable runs
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-        sim.build()
-        variable_names = ["Volume-averaged cell temperature"]
-        out = sim.built_model.export_casadi_objects(variable_names)
+        # Test fails if order not specified
+        with self.assertRaisesRegex(
+            ValueError, "input_parameter_order must be specified"
+        ):
+            model.export_casadi_objects(["a+b"])
+
+        # Fine if order is not specified if there is only one input parameter
+        model = pybamm.BaseModel()
+        p = pybamm.InputParameter("p")
+        model.rhs = {a: -a}
+        model.algebraic = {b: a - b}
+        model.initial_conditions = {a: 1, b: 1}
+        model.variables = {"a+b": a + b - p}
+
+        out = model.export_casadi_objects(["a+b"])
+
+        # Try making a function from the outputs
+        t, x, z, p = out["t"], out["x"], out["z"], out["inputs"]
+        var = out["variables"]["a+b"]
+        var_fn = casadi.Function("var", [t, x, z, p], [var])
+
+        # Test that function values are as expected
+        # a + b - p = 3 + 2 - 7 = -2
+        self.assertEqual(var_fn(6, 3, 2, [7]), -2)
 
         # Test fails if not discretised
+        model = pybamm.lithium_ion.SPMe()
         with self.assertRaisesRegex(
             pybamm.DiscretisationError, "Cannot automatically discretise model"
         ):
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index b2fa082f7e..6ef7a02760 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -44,7 +44,6 @@
 'total interfacial current density as a state': 'false' (possible: ['false', 'true'])
 'working electrode': 'both' (possible: ['both', 'negative', 'positive'])
 'x-average side reactions': 'false' (possible: ['false', 'true'])
-'external submodels': []
 'timescale': 'default'
 """  # noqa: E501
 
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
index af208bd432..1d3a7dd954 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spme.py
@@ -10,22 +10,23 @@ class TestSPMe(BaseUnitTestLithiumIon, unittest.TestCase):
     def setUp(self):
         self.model = pybamm.lithium_ion.SPMe
 
-    def test_external_variables(self):
-        # a concatenation
-        model_options = {"external submodels": ["electrolyte diffusion"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        self.assertEqual(
-            model.external_variables[0],
-            model.variables["Porosity times concentration"],
-        )
-
-        # a variable
-        model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
-        model = pybamm.lithium_ion.SPMe(model_options)
-        self.assertEqual(
-            model.external_variables[0],
-            model.variables["Volume-averaged cell temperature"],
-        )
+    # def test_external_variables(self):
+    #     # To Do: replace external variable with input
+    #     # a concatenation
+    #     model_options = {"external submodels": ["electrolyte diffusion"]}
+    #     model = pybamm.lithium_ion.SPMe(model_options)
+    #     self.assertEqual(
+    #         model.external_variables[0],
+    #         model.variables["Porosity times concentration"],
+    #     )
+
+    #     # a variable
+    #     model_options = {"thermal": "lumped", "external submodels": ["thermal"]}
+    #     model = pybamm.lithium_ion.SPMe(model_options)
+    #     self.assertEqual(
+    #         model.external_variables[0],
+    #         model.variables["Volume-averaged cell temperature"],
+    #     )
 
     def test_electrolyte_options(self):
         options = {"electrolyte conductivity": "full"}
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index 231832e74f..dcb71e51a4 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -149,34 +149,6 @@ def test_set_crate(self):
         self.assertEqual(sim.parameter_values["Current function [A]"], 2 * current_1C)
         self.assertEqual(sim.C_rate, 2)
 
-    def test_set_external_variable(self):
-        model_options = {
-            "thermal": "lumped",
-            "external submodels": ["thermal", "negative primary particle"],
-        }
-        model = pybamm.lithium_ion.SPMe(model_options)
-        sim = pybamm.Simulation(model)
-
-        Nr = model.default_var_pts["r_n"]
-
-        T_av = 0
-        c_s_n_av = np.ones((Nr, 1)) * 0.5
-        external_variables = {
-            "Volume-averaged cell temperature": T_av,
-            "X-averaged negative particle concentration": c_s_n_av,
-        }
-
-        # Step
-        dt = 0.1
-        for _ in range(5):
-            sim.step(dt, external_variables=external_variables)
-        sim.plot(testing=True)
-
-        # Solve
-        t_eval = np.linspace(0, 3600)
-        sim.solve(t_eval, external_variables=external_variables)
-        sim.plot(testing=True)
-
     def test_step(self):
 
         dt = 0.001
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 1174da8f90..439e0b603e 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -443,29 +443,6 @@ def test_model_solver_dae_inputs_in_initial_conditions(self):
             solution.y.full()[-1], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
-    def test_model_solver_with_external(self):
-        # Create model
-        model = pybamm.BaseModel()
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=domain)
-        var2 = pybamm.Variable("var2", domain=domain)
-        model.rhs = {var1: -var2}
-        model.initial_conditions = {var1: 1}
-        model.external_variables = [var2]
-        model.variables = {"var1": var1, "var2": var2}
-        # 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.CasadiSolver(rtol=1e-8, atol=1e-8)
-        t_eval = np.linspace(0, 10, 100)
-        solution = solver.solve(model, t_eval, external_variables={"var2": 0.5})
-        np.testing.assert_allclose(
-            solution.y.full()[0], 1 - 0.5 * solution.t, rtol=1e-06
-        )
-
     def test_model_solver_with_non_identity_mass(self):
         model = pybamm.BaseModel()
         var1 = pybamm.Variable("var1", domain="negative electrode")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 9d186a9072..d483d4bcea 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -702,30 +702,6 @@ def test_model_solver_dae_inputs_in_initial_conditions(self):
             solution.y[-1], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
-    def test_model_solver_dae_with_external(self):
-        # Create model
-        model = pybamm.BaseModel()
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=domain)
-        var2 = pybamm.Variable("var2", domain=domain)
-        model.rhs = {var1: -var2}
-        model.initial_conditions = {var1: 1}
-        model.external_variables = [var2]
-        model.variables = {"var1": var1, "var2": var2}
-        # 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.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
-        t_eval = np.linspace(0, 10, 100)
-        solution = solver.solve(model, t_eval, external_variables={"var2": 0.5})
-        np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06)
-
     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 e3ebeff99e..6358c063d4 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -368,32 +368,6 @@ def test_model_solver_multiple_inputs_jax_format_error(self):
         ):
             solver.solve(model, t_eval, inputs=inputs_list, nproc=2)
 
-    def test_model_solver_with_external(self):
-        # Create model
-        model = pybamm.BaseModel()
-        domain = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=domain)
-        var2 = pybamm.Variable("var2", domain=domain)
-        model.rhs = {var1: -var2}
-        model.initial_conditions = {var1: 1}
-        model.external_variables = [var2]
-        model.variables = {"var2": var2}
-        # 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)
-        t_eval = np.linspace(0, 10, 100)
-        solution = solver.solve(
-            model, t_eval, external_variables={"var2": 0.5 * np.ones(100)}
-        )
-        np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06)
-
     def test_model_solver_with_event_with_casadi(self):
         # Create model
         model = pybamm.BaseModel()