diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 14ebffa96..b6ff10ea6 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -86,6 +86,7 @@ class Discretisation(object):
         self.set_variable_slices(variables)
 
         # set boundary conditions (only need key ids for boundary_conditions)
+        model = self.set_internal_boundary_conditions(model)
         self.bcs = self.process_boundary_conditions(model)
 
         # set up inplace vs not inplace
@@ -97,6 +98,8 @@ class Discretisation(object):
             # create a blank model so that original model is unchanged
             model_disc = pybamm.BaseModel()
 
+        model_disc.bcs = self.bcs
+
         # Process initial condtions
         ics, concat_ics = self.process_initial_conditions(model)
         model_disc.initial_conditions = ics
@@ -167,6 +170,68 @@ class Discretisation(object):
         # reset discretised_symbols
         self._discretised_symbols = {}
 
+    def set_internal_boundary_conditions(self, model):
+        """
+        A method to set the internal boundary conditions for the submodel.
+        These are required to properly calculate the gradient.
+        Note: this method modifies the state of self.boundary_conditions.
+        """
+
+        def boundary_gradient(left_orphan, right_orphan):
+
+            left_domain = left_orphan.domain[0]
+            right_domain = right_orphan.domain[0]
+
+            left_submesh = self._spatial_methods[left_domain].mesh[left_domain][0]
+            right_submesh = self._spatial_methods[right_domain].mesh[right_domain][0]
+
+            x_left = left_submesh.nodes[-1]
+            x_right = right_submesh.nodes[0]
+
+            dy = pybamm.Index(right_orphan, 0) - pybamm.Index(left_orphan, -1)
+            dx = x_right - x_left
+
+            return dy / dx
+
+        bc_key_ids = [key.id for key in list(model.boundary_conditions.keys())]
+
+        internal_bcs = {}
+        for var in model.boundary_conditions.keys():
+            if isinstance(var, pybamm.Concatenation):
+                children = var.children
+
+                first_child = children[0]
+                first_orphan = first_child.new_copy()
+                next_child = children[1]
+                next_orphan = next_child.new_copy()
+
+                lbc = model.boundary_conditions[var]["left"]
+                rbc = (boundary_gradient(first_orphan, next_orphan), "Neumann")
+
+                if first_child.id not in bc_key_ids:
+                    internal_bcs.update({first_child: {"left": lbc, "right": rbc}})
+
+                for i, _ in enumerate(children[1:-1]):
+                    current_child = next_child
+                    current_orphan = next_orphan
+                    next_child = children[i + 2]
+                    next_orphan = next_child.new_copy()
+
+                    lbc = rbc
+                    rbc = (boundary_gradient(current_orphan, next_orphan), "Neumann")
+                    if current_child.id not in bc_key_ids:
+                        internal_bcs.update(
+                            {current_child: {"left": lbc, "right": rbc}}
+                        )
+
+                lbc = rbc
+                rbc = model.boundary_conditions[var]["right"]
+                if children[-1].id not in bc_key_ids:
+                    internal_bcs.update({children[-1]: {"left": lbc, "right": rbc}})
+
+        model.boundary_conditions.update(internal_bcs)
+        return model
+
     def process_initial_conditions(self, model):
         """Discretise model initial_conditions.
 
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 97e383213..ef4536c24 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -144,7 +144,11 @@ class Index(UnaryOperator):
 
     def __init__(self, child, index, name=None):
         self.index = index
-        if isinstance(index, int):
+        if index == -1:
+            self.slice = slice(index, None)
+            if name is None:
+                name = "Index[-1]"
+        elif isinstance(index, int):
             self.slice = slice(index, index + 1)
             if name is None:
                 name = "Index[" + str(index) + "]"
@@ -158,11 +162,17 @@ class Index(UnaryOperator):
         else:
             raise TypeError("index must be integer or slice")
 
-        if self.slice.stop > child.size:
+        if self.slice in (slice(0, 1), slice(-1, None)):
+            pass
+        elif self.slice.stop > child.size:
             raise ValueError("slice size exceeds child size")
 
         super().__init__(name, child)
 
+        # no domain for integer value
+        if isinstance(index, int):
+            self.domain = []
+
     def jac(self, variable):
         """ See :meth:`pybamm.Symbol.jac()`. """
 
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index b02b90924..d3aa235ea 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -273,6 +273,8 @@ class ParameterValues(dict):
             elif isinstance(symbol, pybamm.BoundaryOperator):
                 # BoundaryValue or BoundaryFlux
                 new_symbol = symbol.__class__(new_child, symbol.side)
+            elif isinstance(symbol, pybamm.Index):
+                new_symbol = symbol.__class__(new_child, symbol.index)
             else:
                 new_symbol = symbol.__class__(new_child)
             # ensure domain remains the same

Add handling gradients of subdomain variables
Currently if we try to take the gradient of `phi_e_n` we get a symbol of the wrong size because `phi_e_n` is not in the boundary conditions (and shouldn't be for checking purposes) instead `phi_e` is in the boundary conditions. 
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index 736b991c1..fcd9218ac 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -50,6 +50,75 @@ class TestDFN(unittest.TestCase):
         modeltest = tests.StandardModelTest(model, var_pts=var_pts)
         modeltest.test_all()
 
+    def test_gradient_of_subdomain_var(self):
+        # a hacky test that will be replaced by a more general one in
+        # standard output tests when #481 is completed
+
+        options = {"thermal": None}
+        model = pybamm.lithium_ion.DFN(options)
+        c_e = model.variables["Electrolyte concentration"]
+        c_e_n, c_e_s, c_e_p = c_e.orphans
+        phi_e = model.variables["Electrolyte potential"]
+        phi_e_n, phi_e_s, phi_e_p = phi_e.orphans
+        model.variables.update(
+            {
+                "Test grad c_e": pybamm.grad(c_e),
+                "Test grad c_e_n": pybamm.grad(c_e_n),
+                "Test grad c_e_s": pybamm.grad(c_e_s),
+                "Test grad c_e_p": pybamm.grad(c_e_p),
+                "Test grad phi_e": pybamm.grad(phi_e),
+                "Test grad phi_e_n": pybamm.grad(phi_e_n),
+                "Test grad phi_e_s": pybamm.grad(phi_e_s),
+                "Test grad phi_e_p": pybamm.grad(phi_e_p),
+            }
+        )
+        geometry = model.default_geometry
+        param = model.default_parameter_values
+        param.process_model(model)
+        param.process_geometry(geometry)
+        mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
+        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        disc.process_model(model)
+
+        t_eval = np.linspace(0, 2, 100)
+        solution = model.default_solver.solve(model, t_eval)
+        t = solution.t
+        x_n = disc.mesh["negative electrode"][0].nodes
+        x_s = disc.mesh["separator"][0].nodes
+        x_p = disc.mesh["positive electrode"][0].nodes
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        x = disc.mesh.combine_submeshes(*whole_cell)[0].nodes
+
+        variables = pybamm.post_process_variables(
+            model.variables, solution.t, solution.y, mesh=disc.mesh
+        )
+
+        grad_c_e_n = variables["Test grad c_e_n"]
+        grad_c_e_s = variables["Test grad c_e_s"]
+        grad_c_e_p = variables["Test grad c_e_p"]
+        grad_c_e = variables["Test grad c_e"]
+
+        grad_c_e_combined = np.concatenate(
+            (grad_c_e_n(t, x_n), grad_c_e_s(t, x_s), grad_c_e_p(t, x_p)), axis=0
+        )
+
+        np.testing.assert_array_almost_equal(
+            grad_c_e(t, x), grad_c_e_combined, decimal=6
+        )
+
+        grad_phi_e_n = variables["Test grad phi_e_n"]
+        grad_phi_e_s = variables["Test grad phi_e_s"]
+        grad_phi_e_p = variables["Test grad phi_e_p"]
+        grad_phi_e = variables["Test grad phi_e"]
+
+        grad_phi_e_combined = np.concatenate(
+            (grad_phi_e_n(t, x_n), grad_phi_e_s(t, x_s), grad_phi_e_p(t, x_p)), axis=0
+        )
+
+        np.testing.assert_array_almost_equal(
+            grad_phi_e(t, x), grad_phi_e_combined, decimal=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 964a34010..7357c6df8 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -44,6 +44,23 @@ class TestDiscretise(unittest.TestCase):
         disc = pybamm.Discretisation(None, None)
         self.assertEqual(disc._spatial_methods, {})
 
+    def test_add_internal_boundary_conditions(self):
+        model = pybamm.BaseModel()
+        c_e = pybamm.standard_variables.c_e
+        lbc = (pybamm.Scalar(0), "Neumann")
+        rbc = (pybamm.Scalar(0), "Neumann")
+        model.boundary_conditions = {c_e: {"left": lbc, "right": rbc}}
+
+        mesh = get_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.SpatialMethod,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        model = disc.set_internal_boundary_conditions(model)
+
+        for child in c_e.children:
+            self.assertTrue(child in model.boundary_conditions.keys())
+
     def test_discretise_slicing(self):
         # create discretisation
         mesh = get_mesh_for_testing()

pip install -e .[dev,docs]
parse_log_pytest
pytest
pytest pytest-cov pytest-xdist pytest-mock pytest-asyncio
apt-get update apt-get install -y gcc gfortran libopenblas-dev liblapack-dev
3.9
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.16
anyio==4.9.0
anytree==2.12.1
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
asttokens==3.0.0
async-lru==2.0.5
attrs==25.3.0
autograd==1.7.0
babel==2.17.0
beautifulsoup4==4.13.3
bleach==6.2.0
certifi==2025.1.31
cffi==1.17.1
charset-normalizer==3.4.1
comm==0.2.2
contourpy==1.3.0
coverage==7.8.0
cycler==0.12.1
debugpy==1.8.13
decorator==5.2.1
defusedxml==0.7.1
docutils==0.21.2
exceptiongroup @ file:///croot/exceptiongroup_1706031385326/work
execnet==2.1.1
executing==2.2.0
fastjsonschema==2.21.1
flake8==7.2.0
fonttools==4.56.0
fqdn==1.5.1
guzzle_sphinx_theme==0.7.11
h11==0.14.0
httpcore==1.0.7
httpx==0.28.1
idna==3.10
imagesize==1.4.1
importlib_metadata==8.6.1
importlib_resources==6.5.2
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
ipykernel==6.29.5
ipython==8.18.1
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.10.0
jsonpointer==3.0.0
jsonschema==4.23.0
jsonschema-specifications==2024.10.1
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.12.0
jupyter-lsp==2.2.5
jupyter_client==8.6.3
jupyter_core==5.7.2
jupyter_server==2.15.0
jupyter_server_terminals==0.5.3
jupyterlab==4.3.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.27.3
jupyterlab_widgets==3.0.13
kiwisolver==1.4.7
MarkupSafe==3.0.2
matplotlib==3.9.4
matplotlib-inline==0.1.7
mccabe==0.7.0
mistune==3.1.3
nbclient==0.10.2
nbconvert==7.16.6
nbformat==5.10.4
nest-asyncio==1.6.0
notebook==7.3.3
notebook_shim==0.2.4
numpy==2.0.2
overrides==7.7.0
packaging @ file:///croot/packaging_1734472117206/work
pandas==2.2.3
pandocfilters==1.5.1
parso==0.8.4
pexpect==4.9.0
pillow==11.1.0
platformdirs==4.3.7
pluggy @ file:///croot/pluggy_1733169602837/work
prometheus_client==0.21.1
prompt_toolkit==3.0.50
psutil==7.0.0
ptyprocess==0.7.0
pure_eval==0.2.3
-e git+https://github.com/pybamm-team/PyBaMM.git@3da929593835debf3bb6147d5b68b05a7b13c602#egg=pybamm
pycodestyle==2.13.0
pycparser==2.22
pyflakes==3.3.1
Pygments==2.19.1
pyparsing==3.2.3
pytest @ file:///croot/pytest_1738938843180/work
pytest-asyncio==0.26.0
pytest-cov==6.0.0
pytest-mock==3.14.0
pytest-xdist==3.6.1
python-dateutil==2.9.0.post0
python-json-logger==3.3.0
pytz==2025.2
PyYAML==6.0.2
pyzmq==26.3.0
referencing==0.36.2
requests==2.32.3
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.24.0
scipy==1.13.1
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.6
Sphinx==7.4.7
sphinxcontrib-applehelp==2.0.0
sphinxcontrib-devhelp==2.0.0
sphinxcontrib-htmlhelp==2.1.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==2.0.0
sphinxcontrib-serializinghtml==2.0.0
stack-data==0.6.3
terminado==0.18.1
tinycss2==1.4.0
tomli @ file:///opt/conda/conda-bld/tomli_1657175507142/work
tornado==6.4.2
traitlets==5.14.3
types-python-dateutil==2.9.0.20241206
typing_extensions==4.13.0
tzdata==2025.2
uri-template==1.3.0
urllib3==2.3.0
wcwidth==0.2.13
webcolors==24.11.1
webencodings==0.5.1
websocket-client==1.8.0
widgetsnbextension==4.0.13
zipp==3.21.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - exceptiongroup=1.2.0=py39h06a4308_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=3.0.16=h5eee18b_0
  - packaging=24.2=py39h06a4308_0
  - pip=25.0=py39h06a4308_0
  - pluggy=1.5.0=py39h06a4308_0
  - pytest=8.3.4=py39h06a4308_0
  - python=3.9.21=he870216_1
  - readline=8.2=h5eee18b_0
  - setuptools=75.8.0=py39h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - tomli=2.0.1=py39h06a4308_0
  - wheel=0.45.1=py39h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.16
      - anyio==4.9.0
      - anytree==2.12.1
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - arrow==1.3.0
      - asttokens==3.0.0
      - async-lru==2.0.5
      - attrs==25.3.0
      - autograd==1.7.0
      - babel==2.17.0
      - beautifulsoup4==4.13.3
      - bleach==6.2.0
      - certifi==2025.1.31
      - cffi==1.17.1
      - charset-normalizer==3.4.1
      - comm==0.2.2
      - contourpy==1.3.0
      - coverage==7.8.0
      - cycler==0.12.1
      - debugpy==1.8.13
      - decorator==5.2.1
      - defusedxml==0.7.1
      - docutils==0.21.2
      - execnet==2.1.1
      - executing==2.2.0
      - fastjsonschema==2.21.1
      - flake8==7.2.0
      - fonttools==4.56.0
      - fqdn==1.5.1
      - guzzle-sphinx-theme==0.7.11
      - h11==0.14.0
      - httpcore==1.0.7
      - httpx==0.28.1
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==8.6.1
      - importlib-resources==6.5.2
      - ipykernel==6.29.5
      - ipython==8.18.1
      - ipywidgets==8.1.5
      - isoduration==20.11.0
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.10.0
      - jsonpointer==3.0.0
      - jsonschema==4.23.0
      - jsonschema-specifications==2024.10.1
      - jupyter==1.1.1
      - jupyter-client==8.6.3
      - jupyter-console==6.6.3
      - jupyter-core==5.7.2
      - jupyter-events==0.12.0
      - jupyter-lsp==2.2.5
      - jupyter-server==2.15.0
      - jupyter-server-terminals==0.5.3
      - jupyterlab==4.3.6
      - jupyterlab-pygments==0.3.0
      - jupyterlab-server==2.27.3
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.7
      - markupsafe==3.0.2
      - matplotlib==3.9.4
      - matplotlib-inline==0.1.7
      - mccabe==0.7.0
      - mistune==3.1.3
      - nbclient==0.10.2
      - nbconvert==7.16.6
      - nbformat==5.10.4
      - nest-asyncio==1.6.0
      - notebook==7.3.3
      - notebook-shim==0.2.4
      - numpy==2.0.2
      - overrides==7.7.0
      - pandas==2.2.3
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pexpect==4.9.0
      - pillow==11.1.0
      - platformdirs==4.3.7
      - prometheus-client==0.21.1
      - prompt-toolkit==3.0.50
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pure-eval==0.2.3
      - pycodestyle==2.13.0
      - pycparser==2.22
      - pyflakes==3.3.1
      - pygments==2.19.1
      - pyparsing==3.2.3
      - pytest-asyncio==0.26.0
      - pytest-cov==6.0.0
      - pytest-mock==3.14.0
      - pytest-xdist==3.6.1
      - python-dateutil==2.9.0.post0
      - python-json-logger==3.3.0
      - pytz==2025.2
      - pyyaml==6.0.2
      - pyzmq==26.3.0
      - referencing==0.36.2
      - requests==2.32.3
      - rfc3339-validator==0.1.4
      - rfc3986-validator==0.1.1
      - rpds-py==0.24.0
      - scipy==1.13.1
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - snowballstemmer==2.2.0
      - soupsieve==2.6
      - sphinx==7.4.7
      - sphinxcontrib-applehelp==2.0.0
      - sphinxcontrib-devhelp==2.0.0
      - sphinxcontrib-htmlhelp==2.1.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==2.0.0
      - sphinxcontrib-serializinghtml==2.0.0
      - stack-data==0.6.3
      - terminado==0.18.1
      - tinycss2==1.4.0
      - tornado==6.4.2
      - traitlets==5.14.3
      - types-python-dateutil==2.9.0.20241206
      - typing-extensions==4.13.0
      - tzdata==2025.2
      - uri-template==1.3.0
      - urllib3==2.3.0
      - wcwidth==0.2.13
      - webcolors==24.11.1
      - webencodings==0.5.1
      - websocket-client==1.8.0
      - widgetsnbextension==4.0.13
      - zipp==3.21.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_add_internal_boundary_conditions
tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast_2D tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenate_in_order tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_of_scalars tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_slicing tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_spatial_operator tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_exceptions tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_no_mesh tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_outer tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_complex_expression tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_dict tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_concatenation tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_dae tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_not_inplace tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_ode tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_symbol_base tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_variables_dict
pybamm/parameters/parameter_values.py pybamm/discretisations/discretisation.py pybamm/expression_tree/unary_operators.py
pybamm/parameters/parameter_values.py pybamm/discretisations/discretisation.py pybamm/expression_tree/unary_operators.py