diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 80388c682..bab0dea5a 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -526,10 +526,12 @@ class BinaryOperator(pybamm.Symbol):
             right = pybamm.Scalar(right)
         domain = self.get_children_domains(left.domain, right.domain)
         super().__init__(name, children=[left, right], domain=domain)
+        self.left = self.children[0]
+        self.right = self.children[1]
 
     def __str__(self):
         """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "{!s} {} {!s}".format(self.children[0], self.name, self.children[1])
+        return "{!s} {} {!s}".format(self.left, self.name, self.right)
 
     def get_children_domains(self, ldomain, rdomain):
         if ldomain == rdomain:
@@ -550,8 +552,8 @@ class BinaryOperator(pybamm.Symbol):
 
     def simplify(self):
         """ See :meth:`pybamm.Symbol.simplify()`. """
-        left = self.children[0].simplify()
-        right = self.children[1].simplify()
+        left = self.left.simplify()
+        right = self.right.simplify()
 
         # _binary_simplify defined in derived classes for specific rules
         new_node = self._binary_simplify(left, right)
@@ -569,14 +571,14 @@ class BinaryOperator(pybamm.Symbol):
             try:
                 return known_evals[id], known_evals
             except KeyError:
-                left, known_evals = self.children[0].evaluate(t, y, known_evals)
-                right, known_evals = self.children[1].evaluate(t, y, known_evals)
+                left, known_evals = self.left.evaluate(t, y, known_evals)
+                right, known_evals = self.right.evaluate(t, y, known_evals)
                 value = self._binary_evaluate(left, right)
                 known_evals[id] = value
                 return value, known_evals
         else:
-            left = self.children[0].evaluate(t, y)
-            right = self.children[1].evaluate(t, y)
+            left = self.left.evaluate(t, y)
+            right = self.right.evaluate(t, y)
             return self._binary_evaluate(left, right)
 
 
@@ -656,14 +658,14 @@ class Addition(BinaryOperator):
         if variable.id == self.id:
             return pybamm.Scalar(1)
         else:
-            return self.children[0].diff(variable) + self.children[1].diff(variable)
+            return self.left.diff(variable) + self.right.diff(variable)
 
     def jac(self, variable):
         """ See :meth:`pybamm.Symbol.jac()`. """
         if variable.id == self.id:
             return pybamm.Scalar(1)
         else:
-            return self.children[0].jac(variable) + self.children[1].jac(variable)
+            return self.left.jac(variable) + self.right.jac(variable)
 
     def _binary_evaluate(self, left, right):
         """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
@@ -697,11 +699,11 @@ class Subtraction(BinaryOperator):
         if variable.id == self.id:
             return pybamm.Scalar(1)
         else:
-            return self.children[0].diff(variable) - self.children[1].diff(variable)
+            return self.left.diff(variable) - self.right.diff(variable)
 
     def jac(self, variable):
         """ See :meth:`pybamm.Symbol.jac()`. """
-        return self.children[0].jac(variable) - self.children[1].jac(variable)
+        return self.left.jac(variable) - self.right.jac(variable)
 
     def _binary_evaluate(self, left, right):
         """ See :meth:`pybamm.BinaryOperator._binary_evaluate()`. """
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index c9d5ad7e1..d0af9a374 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -53,16 +53,17 @@ class Concatenation(pybamm.Symbol):
 
     def evaluate(self, t=None, y=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
+        children = self.cached_children
         if known_evals is not None:
             if self.id not in known_evals:
-                children_eval = [None] * len(self.children)
-                for idx, child in enumerate(self.children):
+                children_eval = [None] * len(children)
+                for idx, child in enumerate(children):
                     children_eval[idx], known_evals = child.evaluate(t, y, known_evals)
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
-            children_eval = [None] * len(self.children)
-            for idx, child in enumerate(self.children):
+            children_eval = [None] * len(children)
+            for idx, child in enumerate(children):
                 children_eval[idx] = child.evaluate(t, y)
             return self._concatenation_evaluate(children_eval)
 
@@ -102,14 +103,15 @@ class NumpyConcatenation(Concatenation):
 
     def jac(self, variable):
         """ See :meth:`pybamm.Symbol.jac()`. """
-        if len(self.children) == 0:
+        children = self.cached_children
+        if len(children) == 0:
             return pybamm.Scalar(0)
         else:
-            return SparseStack(*[child.jac(variable) for child in self.children])
+            return SparseStack(*[child.jac(variable) for child in children])
 
     def simplify(self):
         """ See :meth:`pybamm.Symbol.simplify()`. """
-        children = [child.simplify() for child in self.children]
+        children = [child.simplify() for child in self.cached_children]
 
         new_node = self.__class__(*children)
 
@@ -171,7 +173,7 @@ class DomainConcatenation(Concatenation):
 
             # create disc of domain => slice for each child
             self._children_slices = [
-                self.create_slices(child) for child in self.children
+                self.create_slices(child) for child in self.cached_children
             ]
         else:
             self._mesh = copy.copy(copy_this._mesh)
@@ -217,14 +219,15 @@ class DomainConcatenation(Concatenation):
 
     def jac(self, variable):
         """ See :meth:`pybamm.Symbol.jac()`. """
-        if len(self.children) == 0:
+        children = self.cached_children
+        if len(children) == 0:
             return pybamm.Scalar(0)
         else:
-            return SparseStack(*[child.jac(variable) for child in self.children])
+            return SparseStack(*[child.jac(variable) for child in children])
 
     def simplify(self):
         """ See :meth:`pybamm.Symbol.simplify()`. """
-        children = self.children
+        children = self.cached_children
         # Simplify Concatenation of StateVectors to a single StateVector
         if all([isinstance(x, pybamm.StateVector) for x in children]) and all(
             [
@@ -270,16 +273,17 @@ class SparseStack(Concatenation):
 
     def evaluate(self, t=None, y=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
+        children = self.cached_children
         if known_evals is not None:
             if self.id not in known_evals:
-                children_eval = [None] * len(self.children)
-                for idx, child in enumerate(self.children):
+                children_eval = [None] * len(children)
+                for idx, child in enumerate(children):
                     children_eval[idx], known_evals = child.evaluate(t, y, known_evals)
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
-            children_eval = [None] * len(self.children)
-            for idx, child in enumerate(self.children):
+            children_eval = [None] * len(children)
+            for idx, child in enumerate(children):
                 children_eval[idx] = child.evaluate(t, y)
             return self._concatenation_evaluate(children_eval)
 
@@ -292,7 +296,7 @@ class SparseStack(Concatenation):
 
     def simplify(self):
         """ See :meth:`pybamm.Symbol.simplify()`. """
-        children = [child.simplify() for child in self.children]
+        children = [child.simplify() for child in self.cached_children]
 
         new_node = self.__class__(*children)
 
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 7611b986f..978a7d802 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -57,6 +57,9 @@ class Symbol(anytree.NodeMixin):
             # this also adds copy.copy(child) to self.children
             copy.copy(child).parent = self
 
+        # cache children
+        self.cached_children = super(Symbol, self).children
+
         # Set domain (and hence id)
         self.domain = domain
 
@@ -64,6 +67,10 @@ class Symbol(anytree.NodeMixin):
         self._has_left_ghost_cell = False
         self._has_right_ghost_cell = False
 
+    @property
+    def children(self):
+        return self.cached_children
+
     @property
     def name(self):
         """name of the node"""
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 6711fb1cf..52c74fca5 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -31,14 +31,15 @@ class UnaryOperator(pybamm.Symbol):
 
     def __init__(self, name, child):
         super().__init__(name, children=[child], domain=child.domain)
+        self.child = self.children[0]
 
     def __str__(self):
         """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "{}({!s})".format(self.name, self.children[0])
+        return "{}({!s})".format(self.name, self.child)
 
     def simplify(self):
         """ See :meth:`pybamm.Symbol.simplify()`. """
-        child = self.children[0].simplify()
+        child = self.child.simplify()
 
         # _binary_simplify defined in derived classes for specific rules
         new_node = self._unary_simplify(child)
@@ -58,11 +59,11 @@ class UnaryOperator(pybamm.Symbol):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
-                child, known_evals = self.children[0].evaluate(t, y, known_evals)
+                child, known_evals = self.child.evaluate(t, y, known_evals)
                 known_evals[self.id] = self._unary_evaluate(child)
             return known_evals[self.id], known_evals
         else:
-            child = self.children[0].evaluate(t, y)
+            child = self.child.evaluate(t, y)
             return self._unary_evaluate(child)
 
 
@@ -78,18 +79,18 @@ class Negate(UnaryOperator):
 
     def __str__(self):
         """ See :meth:`pybamm.Symbol.__str__()`. """
-        return "{}{!s}".format(self.name, self.children[0])
+        return "{}{!s}".format(self.name, self.child)
 
     def diff(self, variable):
         """ See :meth:`pybamm.Symbol.diff()`. """
         if variable.id == self.id:
             return pybamm.Scalar(1)
         else:
-            return -self.children[0].diff(variable)
+            return -self.child.diff(variable)
 
     def jac(self, variable):
         """ See :meth:`pybamm.Symbol.jac()`. """
-        return -self.children[0].jac(variable)
+        return -self.child.jac(variable)
 
     def _unary_evaluate(self, child):
         """ See :meth:`UnaryOperator._unary_evaluate()`. """
@@ -192,7 +193,7 @@ class Function(UnaryOperator):
             # If self.func() takes no parameters then we can always simplify it
             return pybamm.Scalar(self.func())
         else:
-            child = self.children[0].simplify()
+            child = self.child.simplify()
 
             new_node = pybamm.Function(self.func, child)
 

Cache children
Predefine children, such as `left` and `right` in `BinaryOperator`, to get around `self.children` being slow
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index 40e157886..a929e6698 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -17,6 +17,28 @@ class TestSymbol(unittest.TestCase):
         self.assertEqual(sym.name, "a symbol")
         self.assertEqual(str(sym), "a symbol")
 
+    def test_cached_children(self):
+        symc1 = pybamm.Symbol("child1")
+        symc2 = pybamm.Symbol("child2")
+        symc3 = pybamm.Symbol("child3")
+        symp = pybamm.Symbol("parent", children=[symc1, symc2])
+
+        # test tuples of children for equality based on their name
+        def check_are_equal(children1, children2):
+            self.assertEqual(len(children1), len(children2))
+            for i in range(len(children1)):
+                self.assertEqual(children1[i].name, children2[i].name)
+
+        check_are_equal(symp.children, super(pybamm.Symbol, symp).children)
+        check_are_equal(symp.children, (symc1, symc2))
+
+        # update children, since we cache the children they will be unchanged
+        symc3.parent = symp
+        check_are_equal(symp.children, (symc1, symc2))
+
+        # check that the *actual* children are updated
+        check_are_equal(super(pybamm.Symbol, symp).children, (symc1, symc2, symc3))
+
     def test_symbol_simplify(self):
         a = pybamm.Scalar(0)
         b = pybamm.Scalar(1)

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest pytest-cov pytest-xdist pytest-mock pytest-asyncio pytest
apt-get update apt-get install -y python3-dev gfortran gcc libopenblas-dev liblapack-dev
3.6
requirements/base.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs==22.2.0
autograd==1.6.2
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
coverage==6.2
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
execnet==1.9.0
flake8==5.0.4
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig==1.1.1
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging==21.3
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy==1.0.0
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py==1.11.0
-e git+https://github.com/pybamm-team/PyBaMM.git@2fa950247014390bb6fd73f8ae9def38ccea97ed#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing==3.1.4
pyrsistent==0.18.0
pytest==7.0.1
pytest-asyncio==0.16.0
pytest-cov==4.0.0
pytest-mock==3.6.1
pytest-xdist==3.0.2
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
tomli==1.2.3
tornado==6.1
traitlets==4.3.3
typing_extensions==4.1.1
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp==3.6.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
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - 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=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - attrs==22.2.0
      - autograd==1.6.2
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - coverage==6.2
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - execnet==1.9.0
      - flake8==5.0.4
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - iniconfig==1.1.1
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - pluggy==1.0.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - py==1.11.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyparsing==3.1.4
      - pyrsistent==0.18.0
      - pytest==7.0.1
      - pytest-asyncio==0.16.0
      - pytest-cov==4.0.0
      - pytest-mock==3.6.1
      - pytest-xdist==3.0.2
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tomli==1.2.3
      - tornado==6.1
      - traitlets==4.3.3
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
      - zipp==3.6.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_cached_children
tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_visualise
tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_ghost_cell_flaf tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_has_spatial_derivatives tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_multiple_symbols tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_orphans tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_diff tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_domains tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluates_to_number tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluation tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_init tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_is_constant tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_methods tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_repr tests/unit/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_simplify
pybamm/expression_tree/unary_operators.py pybamm/expression_tree/binary_operators.py pybamm/expression_tree/symbol.py pybamm/expression_tree/concatenations.py
pybamm/expression_tree/unary_operators.py pybamm/expression_tree/binary_operators.py pybamm/expression_tree/symbol.py pybamm/expression_tree/concatenations.py