diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0218499aa..58ce82fd8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added `BackwardIndefiniteIntegral` symbol ([#1014](https://github.com/pybamm-team/PyBaMM/pull/1014))
 -   Added SEI film resistance as an option ([#994](https://github.com/pybamm-team/PyBaMM/pull/994))
 -   Added tab, edge, and surface cooling ([#965](https://github.com/pybamm-team/PyBaMM/pull/965))
 -   Added functionality to solver to automatically discretise a 0D model ([#947](https://github.com/pybamm-team/PyBaMM/pull/947))
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index e9cae42dd..0d6d2f622 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -836,7 +836,13 @@ class Discretisation(object):
                 return child_spatial_method.boundary_mass_matrix(child, self.bcs)
 
             elif isinstance(symbol, pybamm.IndefiniteIntegral):
-                return child_spatial_method.indefinite_integral(child, disc_child)
+                return child_spatial_method.indefinite_integral(
+                    child, disc_child, "forward"
+                )
+            elif isinstance(symbol, pybamm.BackwardIndefiniteIntegral):
+                return child_spatial_method.indefinite_integral(
+                    child, disc_child, "backward"
+                )
 
             elif isinstance(symbol, pybamm.Integral):
                 out = child_spatial_method.integral(child, disc_child)
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index b1d785c39..78a08181a 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -510,14 +510,8 @@ class Integral(SpatialOperator):
         return False
 
 
-class IndefiniteIntegral(Integral):
-    """A node in the expression tree representing an indefinite integral operator
-
-    .. math::
-        I = \\int_{x_\text{min}}^{x}\\!f(u)\\,du
-
-    where :math:`u\\in\\text{domain}` which can represent either a
-    spatial or temporal variable.
+class BaseIndefiniteIntegral(Integral):
+    """Base class for indefinite integrals (forward or backward).
 
     Parameters
     ----------
@@ -540,15 +534,73 @@ class IndefiniteIntegral(Integral):
         super().__init__(child, integration_variable)
         # overwrite domains with child domains
         self.copy_domains(child)
+
+    def _evaluate_for_shape(self):
+        return self.children[0].evaluate_for_shape()
+
+    def evaluates_on_edges(self):
+        # If child evaluates on edges, indefinite integral doesn't
+        # If child doesn't evaluate on edges, indefinite integral does
+        return not self.child.evaluates_on_edges()
+
+
+class IndefiniteIntegral(BaseIndefiniteIntegral):
+    """A node in the expression tree representing an indefinite integral operator
+
+    .. math::
+        I = \\int_{x_\text{min}}^{x}\\!f(u)\\,du
+
+    where :math:`u\\in\\text{domain}` which can represent either a
+    spatial or temporal variable.
+
+    Parameters
+    ----------
+    function : :class:`pybamm.Symbol`
+        The function to be integrated (will become self.children[0])
+    integration_variable : :class:`pybamm.IndependentVariable`
+        The variable over which to integrate
+
+    **Extends:** :class:`BaseIndefiniteIntegral`
+    """
+
+    def __init__(self, child, integration_variable):
+        super().__init__(child, integration_variable)
         # Overwrite the name
         self.name = "{} integrated w.r.t {}".format(
-            child.name, integration_variable.name
+            child.name, self.integration_variable[0].name
         )
         if isinstance(integration_variable, pybamm.SpatialVariable):
-            self.name += " on {}".format(integration_variable.domain)
+            self.name += " on {}".format(self.integration_variable[0].domain)
 
-    def _evaluate_for_shape(self):
-        return self.children[0].evaluate_for_shape()
+
+class BackwardIndefiniteIntegral(BaseIndefiniteIntegral):
+    """A node in the expression tree representing a backward indefinite integral
+    operator
+
+    .. math::
+        I = \\int_{x}^{x_\text{max}}\\!f(u)\\,du
+
+    where :math:`u\\in\\text{domain}` which can represent either a
+    spatial or temporal variable.
+
+    Parameters
+    ----------
+    function : :class:`pybamm.Symbol`
+        The function to be integrated (will become self.children[0])
+    integration_variable : :class:`pybamm.IndependentVariable`
+        The variable over which to integrate
+
+    **Extends:** :class:`BaseIndefiniteIntegral`
+    """
+
+    def __init__(self, child, integration_variable):
+        super().__init__(child, integration_variable)
+        # Overwrite the name
+        self.name = "{} integrated backward w.r.t {}".format(
+            child.name, self.integration_variable[0].name
+        )
+        if isinstance(integration_variable, pybamm.SpatialVariable):
+            self.name += " on {}".format(self.integration_variable[0].domain)
 
 
 class DefiniteIntegralVector(SpatialOperator):
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index d748f5590..671b388c0 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -5,6 +5,7 @@ import pybamm
 
 from scipy.sparse import (
     diags,
+    spdiags,
     eye,
     kron,
     csr_matrix,
@@ -248,6 +249,8 @@ class FiniteVolume(pybamm.SpatialMethod):
         else:
             out = integration_vector @ discretised_child
 
+        out.copy_domains(child)
+
         return out
 
     def definite_integral_matrix(self, domain, vector_type="row"):
@@ -296,15 +299,19 @@ class FiniteVolume(pybamm.SpatialMethod):
         matrix = csr_matrix(kron(eye(second_dim_len), vector))
         return pybamm.Matrix(matrix)
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """Implementation of the indefinite integral operator. """
 
         # Different integral matrix depending on whether the integrand evaluates on
         # edges or nodes
         if child.evaluates_on_edges():
-            integration_matrix = self.indefinite_integral_matrix_edges(child.domain)
+            integration_matrix = self.indefinite_integral_matrix_edges(
+                child.domain, direction
+            )
         else:
-            integration_matrix = self.indefinite_integral_matrix_nodes(child.domain)
+            integration_matrix = self.indefinite_integral_matrix_nodes(
+                child.domain, direction
+            )
 
         # Don't need to check for spherical domains as spherical polars
         # only change the diveregence (childs here have grad and no div)
@@ -314,10 +321,28 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return out
 
-    def indefinite_integral_matrix_edges(self, domain):
+    def indefinite_integral_matrix_edges(self, domain, direction):
         """
         Matrix for finite-volume implementation of the indefinite integral where the
-        integrand is evaluated on mesh edges
+        integrand is evaluated on mesh edges (shape (n+1, 1)).
+        The integral will then be evaluated on mesh nodes (shape (n, 1)).
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) of integration
+        direction : str
+            The direction of integration (forward or backward). See notes.
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The finite volume integral matrix for the domain
+
+        Notes
+        -----
+
+        **Forward integral**
 
         .. math::
             F(x) = \\int_0^x\\!f(u)\\,du
@@ -335,16 +360,81 @@ class FiniteVolume(pybamm.SpatialMethod):
         Hence we must have
 
         - :math:`F_0 = du_{1/2} * f_{1/2} / 2`
-        - :math:`F_{i+1} = F_i + du * f_{i+1/2}`
+        - :math:`F_{i+1} = F_i + du_{i+1/2} * f_{i+1/2}`
+
+        Note that :math:`f_{-1/2}` and :math:`f_{end+1/2}` are included in the discrete
+        integrand vector `f`, so we add a column of zeros at each end of the
+        indefinite integral matrix to ignore these.
+
+        **Backward integral**
+
+        .. math::
+            F(x) = \\int_x^end\\!f(u)\\,du
 
-        Note that :math:`f_{-1/2}` and :math:`f_{n+1/2}` are included in the discrete
+        The indefinite integral must satisfy the following conditions:
+
+        - :math:`F(end) = 0`
+        - :math:`f(x) = -\\frac{dF}{dx}`
+
+        or, in discrete form,
+
+        - `BoundaryValue(F, "right") = 0`, i.e. :math:`3*F_{end} - F_{end-1} = 0`
+        - :math:`f_{i+1/2} = -(F_{i+1} - F_i) / dx_{i+1/2}`
+
+        Hence we must have
+
+        - :math:`F_{end} = du_{end+1/2} * f_{end-1/2} / 2`
+        - :math:`F_{i-1} = F_i + du_{i-1/2} * f_{i-1/2}`
+
+        Note that :math:`f_{-1/2}` and :math:`f_{end+1/2}` are included in the discrete
         integrand vector `f`, so we add a column of zeros at each end of the
         indefinite integral matrix to ignore these.
+        """
+
+        # Create appropriate submesh by combining submeshes in domain
+        submesh_list = self.mesh.combine_submeshes(*domain)
+        submesh = submesh_list[0]
+        n = submesh.npts
+        sec_pts = len(submesh_list)
+
+        du_n = submesh.d_nodes
+        if direction == "forward":
+            du_entries = [du_n] * (n - 1)
+            offset = -np.arange(1, n, 1)
+            main_integral_matrix = spdiags(du_entries, offset, n, n - 1)
+            bc_offset_matrix = lil_matrix((n, n - 1))
+            bc_offset_matrix[:, 0] = du_n[0] / 2
+        elif direction == "backward":
+            du_entries = [du_n] * (n + 1)
+            offset = np.arange(n, -1, -1)
+            main_integral_matrix = spdiags(du_entries, offset, n, n - 1)
+            bc_offset_matrix = lil_matrix((n, n - 1))
+            bc_offset_matrix[:, -1] = du_n[-1] / 2
+        sub_matrix = main_integral_matrix + bc_offset_matrix
+        # add a column of zeros at each end
+        zero_col = csr_matrix((n, 1))
+        sub_matrix = hstack([zero_col, sub_matrix, zero_col])
+        # Convert to csr_matrix so that we can take the index (row-slicing), which is
+        # not supported by the default kron format
+        # Note that this makes column-slicing inefficient, but this should not be an
+        # issue
+        matrix = csr_matrix(kron(eye(sec_pts), sub_matrix))
+
+        return pybamm.Matrix(matrix)
+
+    def indefinite_integral_matrix_nodes(self, domain, direction):
+        """
+        Matrix for finite-volume implementation of the (backward) indefinite integral
+        where the integrand is evaluated on mesh nodes (shape (n, 1)).
+        The integral will then be evaluated on mesh edges (shape (n+1, 1)).
+        This is just a straightforward (backward) cumulative sum of the integrand
 
         Parameters
         ----------
         domain : list
             The domain(s) of integration
+        direction : str
+            The direction of integration (forward or backward)
 
         Returns
         -------
@@ -358,16 +448,13 @@ class FiniteVolume(pybamm.SpatialMethod):
         n = submesh.npts
         sec_pts = len(submesh_list)
 
-        du_n = submesh.d_nodes
-        du_entries = [du_n] * (n - 1)
-        offset = -np.arange(1, n, 1)
-        main_integral_matrix = diags(du_entries, offset, shape=(n, n - 1))
-        bc_offset_matrix = lil_matrix((n, n - 1))
-        bc_offset_matrix[:, 0] = du_n[0] / 2
-        sub_matrix = main_integral_matrix + bc_offset_matrix
-        # add a column of zeros at each end
-        zero_col = csr_matrix((n, 1))
-        sub_matrix = hstack([zero_col, sub_matrix, zero_col])
+        du_n = submesh.d_edges
+        du_entries = [du_n] * n
+        if direction == "forward":
+            offset = -np.arange(1, n + 1, 1)  # from -1 down to -n
+        elif direction == "backward":
+            offset = np.arange(n - 1, -1, -1)  # from n-1 down to 0
+        sub_matrix = spdiags(du_entries, offset, n + 1, n)
         # Convert to csr_matrix so that we can take the index (row-slicing), which is
         # not supported by the default kron format
         # Note that this makes column-slicing inefficient, but this should not be an
@@ -473,41 +560,6 @@ class FiniteVolume(pybamm.SpatialMethod):
 
         return dy / dx
 
-    def indefinite_integral_matrix_nodes(self, domain):
-        """
-        Matrix for finite-volume implementation of the indefinite integral where the
-        integrand is evaluated on mesh nodes.
-        This is just a straightforward cumulative sum of the integrand
-
-        Parameters
-        ----------
-        domain : list
-            The domain(s) of integration
-
-        Returns
-        -------
-        :class:`pybamm.Matrix`
-            The finite volume integral matrix for the domain
-        """
-
-        # Create appropriate submesh by combining submeshes in domain
-        submesh_list = self.mesh.combine_submeshes(*domain)
-        submesh = submesh_list[0]
-        n = submesh.npts
-        sec_pts = len(submesh_list)
-
-        du_n = submesh.d_edges
-        du_entries = [du_n] * (n)
-        offset = -np.arange(1, n + 1, 1)
-        sub_matrix = diags(du_entries, offset, shape=(n + 1, n))
-        # Convert to csr_matrix so that we can take the index (row-slicing), which is
-        # not supported by the default kron format
-        # Note that this makes column-slicing inefficient, but this should not be an
-        # issue
-        matrix = csr_matrix(kron(eye(sec_pts), sub_matrix))
-
-        return pybamm.Matrix(matrix)
-
     def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
         """
         Add ghost nodes to a symbol.
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 371d50ebc..41d26539b 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -242,7 +242,7 @@ class SpatialMethod:
         """
         raise NotImplementedError
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """
         Implements the indefinite integral for a spatial method.
 
@@ -252,6 +252,8 @@ class SpatialMethod:
             The symbol to which is being integrated
         discretised_child: :class:`pybamm.Symbol`
             The discretised symbol of the correct size
+        direction : str
+            The direction of integration
 
         Returns
         -------
diff --git a/pybamm/spatial_methods/zero_dimensional_method.py b/pybamm/spatial_methods/zero_dimensional_method.py
index aada27398..845270099 100644
--- a/pybamm/spatial_methods/zero_dimensional_method.py
+++ b/pybamm/spatial_methods/zero_dimensional_method.py
@@ -37,11 +37,16 @@ class ZeroDimensionalSpatialMethod(pybamm.SpatialMethod):
         """
         return pybamm.Matrix(np.ones((1, 1)))
 
-    def indefinite_integral(self, child, discretised_child):
+    def indefinite_integral(self, child, discretised_child, direction):
         """
-        Calculates the zero-dimensional indefinite integral, i.e. the identity operator
+        Calculates the zero-dimensional indefinite integral.
+        If 'direction' is forward, this is the identity operator.
+        If 'direction' is backward, this is the negation operator.
         """
-        return discretised_child
+        if direction == "forward":
+            return discretised_child
+        elif direction == "backward":
+            return -discretised_child
 
     def integral(self, child, discretised_child):
         """

backward indefinite integral
Add new node `BackwardIndefiniteIntegral` for integrating from x to end.

The relationship
```
IndefiniteIntegral(var, x) + BackwardIndefiniteIntegral(var, x) = Integral(var, x)
```
should hold for all x
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 a72faf649..e7ddbf305 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -12,6 +12,7 @@ class TestCopy(unittest.TestCase):
         a = pybamm.Scalar(0)
         b = pybamm.Scalar(1)
         v_n = pybamm.Variable("v", "negative electrode")
+        x_n = pybamm.standard_spatial_vars.x_n
         v_s = pybamm.Variable("v", "separator")
         vec = pybamm.Vector(np.array([1, 2, 3, 4, 5]))
         mesh = get_mesh_for_testing()
@@ -29,6 +30,8 @@ class TestCopy(unittest.TestCase):
             pybamm.grad(v_n),
             pybamm.div(pybamm.grad(v_n)),
             pybamm.Integral(a, pybamm.t),
+            pybamm.IndefiniteIntegral(v_n, x_n),
+            pybamm.BackwardIndefiniteIntegral(v_n, x_n),
             pybamm.BoundaryValue(v_n, "right"),
             pybamm.BoundaryGradient(v_n, "right"),
             pybamm.PrimaryBroadcast(a, "domain"),
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 28d206b92..030deaa73 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -173,6 +173,11 @@ class TestUnaryOperators(unittest.TestCase):
         self.assertEqual(
             inta_sec.auxiliary_domains, {"secondary": ["current collector"]}
         )
+        # backward indefinite integral
+        inta = pybamm.BackwardIndefiniteIntegral(a, x)
+        self.assertEqual(
+            inta.name, "a integrated backward w.r.t x on ['negative electrode']"
+        )
 
         # expected errors
         a = pybamm.Symbol("a", domain=["negative electrode"])
diff --git a/tests/unit/test_spatial_methods/test_base_spatial_method.py b/tests/unit/test_spatial_methods/test_base_spatial_method.py
index 475c3aa91..211b4586c 100644
--- a/tests/unit/test_spatial_methods/test_base_spatial_method.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -24,7 +24,7 @@ class TestSpatialMethod(unittest.TestCase):
         with self.assertRaises(NotImplementedError):
             spatial_method.integral(None, None)
         with self.assertRaises(NotImplementedError):
-            spatial_method.indefinite_integral(None, None)
+            spatial_method.indefinite_integral(None, None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.boundary_integral(None, None, None)
         with self.assertRaises(NotImplementedError):
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 8c261c1f9..c023029b7 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -872,6 +872,55 @@ class TestFiniteVolume(unittest.TestCase):
             left_boundary_value_disc.evaluate(y=c_exact), 0
         )
 
+    def test_backward_indefinite_integral(self):
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # --------------------------------------------------------------------
+        # region which doesn't start at zero
+        phi = pybamm.Variable("phi", domain=["negative electrode", "separator"])
+        i = pybamm.grad(phi)  # create test current (variable on edges)
+        x = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+        int_grad_phi = pybamm.BackwardIndefiniteIntegral(i, x)
+        disc.set_variable_slices([phi])  # i is not a fundamental variable
+        disc._bcs = {
+            phi.id: {
+                "left": (pybamm.Scalar(0), "Neumann"),
+                "right": (pybamm.Scalar(0), "Neumann"),
+            }
+        }
+        int_grad_phi_disc = disc.process_symbol(int_grad_phi)
+        right_boundary_value = pybamm.BoundaryValue(int_grad_phi, "right")
+        right_boundary_value_disc = disc.process_symbol(right_boundary_value)
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+
+        # Test that the backward_integral(grad(phi)) = -phi
+        # constant case
+        phi_exact = np.ones((combined_submesh[0].npts, 1))
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
+        phi_approx += 1  # add constant of integration
+        np.testing.assert_array_equal(phi_exact, phi_approx)
+        self.assertEqual(right_boundary_value_disc.evaluate(y=phi_exact), 0)
+
+        # linear case
+        phi_exact = combined_submesh[0].nodes - combined_submesh[0].edges[-1]
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(phi_exact, -phi_approx)
+        np.testing.assert_array_almost_equal(
+            right_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
+
+        # sine case
+        phi_exact = np.sin(combined_submesh[0].nodes - combined_submesh[0].edges[-1])
+        phi_approx = int_grad_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(phi_exact, -phi_approx)
+        np.testing.assert_array_almost_equal(
+            right_boundary_value_disc.evaluate(y=phi_exact), 0
+        )
+
     def test_indefinite_integral_of_broadcasted_to_cell_edges(self):
         # create discretisation
         mesh = get_mesh_for_testing()
@@ -925,19 +974,86 @@ class TestFiniteVolume(unittest.TestCase):
         # constant case
         phi_exact = np.ones((combined_submesh[0].npts, 1))
         int_phi_exact = combined_submesh[0].edges
-        int_phi_approx = int_phi_disc.evaluate(None, phi_exact)[:, 0]
+        int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_equal(int_phi_exact, int_phi_approx)
         # linear case
-        phi_exact = combined_submesh[0].nodes[:, np.newaxis]
-        int_phi_exact = combined_submesh[0].edges[:, np.newaxis] ** 2 / 2
-        int_phi_approx = int_phi_disc.evaluate(None, phi_exact)
+        phi_exact = combined_submesh[0].nodes
+        int_phi_exact = combined_submesh[0].edges ** 2 / 2
+        int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(int_phi_exact, int_phi_approx)
         # cos case
-        phi_exact = np.cos(combined_submesh[0].nodes[:, np.newaxis])
-        int_phi_exact = np.sin(combined_submesh[0].edges[:, np.newaxis])
-        int_phi_approx = int_phi_disc.evaluate(None, phi_exact)
+        phi_exact = np.cos(combined_submesh[0].nodes)
+        int_phi_exact = np.sin(combined_submesh[0].edges)
+        int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(int_phi_exact, int_phi_approx, decimal=5)
 
+    def test_backward_indefinite_integral_on_nodes(self):
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        phi = pybamm.Variable("phi", domain=["negative electrode", "separator"])
+        x = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
+
+        back_int_phi = pybamm.BackwardIndefiniteIntegral(phi, x)
+        disc.set_variable_slices([phi])
+        back_int_phi_disc = disc.process_symbol(back_int_phi)
+
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        edges = combined_submesh[0].edges
+
+        # constant case
+        phi_exact = np.ones((combined_submesh[0].npts, 1))
+        back_int_phi_exact = edges[-1] - edges
+        back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(back_int_phi_exact, back_int_phi_approx)
+        # linear case
+        phi_exact = combined_submesh[0].nodes
+        back_int_phi_exact = edges[-1] ** 2 / 2 - edges ** 2 / 2
+        back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(back_int_phi_exact, back_int_phi_approx)
+        # cos case
+        phi_exact = np.cos(combined_submesh[0].nodes)
+        back_int_phi_exact = np.sin(edges[-1]) - np.sin(edges)
+        back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
+        np.testing.assert_array_almost_equal(
+            back_int_phi_exact, back_int_phi_approx, decimal=5
+        )
+
+    def test_forward_plus_backward_integral(self):
+        # Test that forward integral + backward integral = integral
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # On nodes
+        phi = pybamm.Variable("phi", domain=["separator", "positive electrode"])
+        x = pybamm.SpatialVariable("x", ["separator", "positive electrode"])
+
+        disc.set_variable_slices([phi])
+
+        full_int_phi = pybamm.PrimaryBroadcastToEdges(
+            pybamm.Integral(phi, x), ["separator", "positive electrode"]
+        )
+        full_int_phi_disc = disc.process_symbol(full_int_phi)
+        int_plus_back_int_phi = pybamm.IndefiniteIntegral(
+            phi, x
+        ) + pybamm.BackwardIndefiniteIntegral(phi, x)
+        int_plus_back_int_phi_disc = disc.process_symbol(int_plus_back_int_phi)
+
+        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+
+        # test
+        for phi_exact in [
+            np.ones((combined_submesh[0].npts, 1)),
+            combined_submesh[0].nodes,
+            np.cos(combined_submesh[0].nodes),
+        ]:
+            np.testing.assert_array_almost_equal(
+                full_int_phi_disc.evaluate(y=phi_exact).flatten(),
+                int_plus_back_int_phi_disc.evaluate(y=phi_exact).flatten(),
+            )
+
     def test_discretise_spatial_variable(self):
         # create discretisation
         mesh = get_mesh_for_testing()
diff --git a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py b/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
index f3ffdd80f..de893d981 100644
--- a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
+++ b/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
@@ -16,8 +16,11 @@ class TestZeroDimensionalSpatialMethod(unittest.TestCase):
 
         a = pybamm.Symbol("a")
         self.assertEqual(a, spatial_method.integral(None, a))
-        self.assertEqual(a, spatial_method.indefinite_integral(None, a))
+        self.assertEqual(a, spatial_method.indefinite_integral(None, a, "forward"))
         self.assertEqual(a, spatial_method.boundary_value_or_flux(None, a))
+        self.assertEqual(
+            (-a).id, spatial_method.indefinite_integral(None, a, "backward").id
+        )
 
         mass_matrix = spatial_method.mass_matrix(None, None)
         self.assertIsInstance(mass_matrix, pybamm.Matrix)

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y git graphviz libopenblas-dev
3.7
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
aiofiles==22.1.0
aiosqlite==0.19.0
alabaster==0.7.13
anyio==3.7.1
anytree==2.12.1
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.2.3
attrs==24.2.0
autograd==1.6.2
Babel==2.14.0
backcall==0.2.0
beautifulsoup4==4.13.3
black==23.3.0
bleach==6.0.0
cached-property==1.5.2
casadi==3.7.0
certifi @ file:///croot/certifi_1671487769961/work/certifi
cffi==1.15.1
charset-normalizer==3.4.1
click==8.1.8
comm==0.1.4
cycler==0.11.0
debugpy==1.7.0
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
exceptiongroup==1.2.2
fastjsonschema==2.21.1
flake8==5.0.4
fonttools==4.38.0
fqdn==1.5.1
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
importlib-metadata==4.2.0
importlib-resources==5.12.0
iniconfig==2.0.0
ipykernel==6.16.2
ipython==7.34.0
ipython-genutils==0.2.0
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.9.16
jsonpointer==3.0.0
jsonschema==4.17.3
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.6.3
jupyter-server==1.24.0
jupyter-ydoc==0.2.5
jupyter_client==7.4.9
jupyter_core==4.12.0
jupyter_server_fileid==0.9.3
jupyter_server_ydoc==0.8.0
jupyterlab==3.5.3
jupyterlab-pygments==0.2.2
jupyterlab_server==2.15.2
jupyterlab_widgets==3.0.13
kiwisolver==1.4.5
MarkupSafe==2.1.5
matplotlib==3.5.3
matplotlib-inline==0.1.6
mccabe==0.7.0
mistune==3.0.2
mypy-extensions==1.0.0
nbclassic==1.2.0
nbclient==0.7.4
nbconvert==7.6.0
nbformat==5.8.0
nest-asyncio==1.6.0
notebook==6.5.7
notebook_shim==0.2.4
numpy==1.21.6
packaging==24.0
pandas==1.3.5
pandocfilters==1.5.1
parso==0.8.4
pathspec==0.11.2
pexpect==4.9.0
pickleshare==0.7.5
Pillow==9.5.0
pkgutil_resolve_name==1.3.10
platformdirs==4.0.0
pluggy==1.2.0
prometheus-client==0.17.1
prompt_toolkit==3.0.48
psutil==7.0.0
ptyprocess==0.7.0
-e git+https://github.com/pybamm-team/PyBaMM.git@a45f180b8d55d00882c4f56d64c557cdf97e819e#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.17.2
pyparsing==3.1.4
pyrsistent==0.19.3
pytest==7.4.4
python-dateutil==2.9.0.post0
python-json-logger==3.0.1
pytz==2025.2
PyYAML==6.0.1
pyzmq==26.2.1
requests==2.31.0
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
scikit-fem==8.1.0
scipy==1.7.3
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.4.1
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.17.1
tinycss2==1.2.1
tomli==2.0.1
tornado==6.2
traitlets==5.9.0
typed-ast==1.5.5
typing_extensions==4.7.1
uri-template==1.3.0
urllib3==2.0.7
wcwidth==0.2.13
webcolors==1.13
webencodings==0.5.1
websocket-client==1.6.1
widgetsnbextension==4.0.13
y-py==0.6.2
ypy-websocket==0.8.4
zipp==3.15.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=2022.12.7=py37h06a4308_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=1.1.1w=h7f8727e_0
  - pip=22.3.1=py37h06a4308_0
  - python=3.7.16=h7a1cb2a_0
  - readline=8.2=h5eee18b_0
  - setuptools=65.6.3=py37h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.38.4=py37h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - aiofiles==22.1.0
      - aiosqlite==0.19.0
      - alabaster==0.7.13
      - anyio==3.7.1
      - anytree==2.12.1
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - arrow==1.2.3
      - attrs==24.2.0
      - autograd==1.6.2
      - babel==2.14.0
      - backcall==0.2.0
      - beautifulsoup4==4.13.3
      - black==23.3.0
      - bleach==6.0.0
      - cached-property==1.5.2
      - casadi==3.7.0
      - cffi==1.15.1
      - charset-normalizer==3.4.1
      - click==8.1.8
      - comm==0.1.4
      - cycler==0.11.0
      - debugpy==1.7.0
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - exceptiongroup==1.2.2
      - fastjsonschema==2.21.1
      - flake8==5.0.4
      - fonttools==4.38.0
      - fqdn==1.5.1
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==4.2.0
      - importlib-resources==5.12.0
      - iniconfig==2.0.0
      - ipykernel==6.16.2
      - ipython==7.34.0
      - ipython-genutils==0.2.0
      - ipywidgets==8.1.5
      - isoduration==20.11.0
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.9.16
      - jsonpointer==3.0.0
      - jsonschema==4.17.3
      - jupyter==1.1.1
      - jupyter-client==7.4.9
      - jupyter-console==6.6.3
      - jupyter-core==4.12.0
      - jupyter-events==0.6.3
      - jupyter-server==1.24.0
      - jupyter-server-fileid==0.9.3
      - jupyter-server-ydoc==0.8.0
      - jupyter-ydoc==0.2.5
      - jupyterlab==3.5.3
      - jupyterlab-pygments==0.2.2
      - jupyterlab-server==2.15.2
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.5
      - markupsafe==2.1.5
      - matplotlib==3.5.3
      - matplotlib-inline==0.1.6
      - mccabe==0.7.0
      - mistune==3.0.2
      - mypy-extensions==1.0.0
      - nbclassic==1.2.0
      - nbclient==0.7.4
      - nbconvert==7.6.0
      - nbformat==5.8.0
      - nest-asyncio==1.6.0
      - notebook==6.5.7
      - notebook-shim==0.2.4
      - numpy==1.21.6
      - packaging==24.0
      - pandas==1.3.5
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pathspec==0.11.2
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==9.5.0
      - pkgutil-resolve-name==1.3.10
      - platformdirs==4.0.0
      - pluggy==1.2.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.48
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.17.2
      - pyparsing==3.1.4
      - pyrsistent==0.19.3
      - pytest==7.4.4
      - python-dateutil==2.9.0.post0
      - python-json-logger==3.0.1
      - pytz==2025.2
      - pyyaml==6.0.1
      - pyzmq==26.2.1
      - requests==2.31.0
      - rfc3339-validator==0.1.4
      - rfc3986-validator==0.1.1
      - scikit-fem==8.1.0
      - scipy==1.7.3
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - snowballstemmer==2.2.0
      - soupsieve==2.4.1
      - 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.17.1
      - tinycss2==1.2.1
      - tomli==2.0.1
      - tornado==6.2
      - traitlets==5.9.0
      - typed-ast==1.5.5
      - typing-extensions==4.7.1
      - uri-template==1.3.0
      - urllib3==2.0.7
      - wcwidth==0.2.13
      - webcolors==1.13
      - webencodings==0.5.1
      - websocket-client==1.6.1
      - widgetsnbextension==4.0.13
      - y-py==0.6.2
      - ypy-websocket==0.8.4
      - zipp==3.15.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_expression_tree/test_operations/test_copy.py::TestCopy::test_symbol_new_copy tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_integral tests/unit/test_spatial_methods/test_base_spatial_method.py::TestSpatialMethod::test_basics tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_backward_indefinite_integral tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_backward_indefinite_integral_on_nodes tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_forward_plus_backward_integral tests/unit/test_spatial_methods/test_zero_dimensional_method.py::TestZeroDimensionalSpatialMethod::test_identity_ops
tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_absolute tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_operators tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_boundary_value tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_delta_function tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_diff tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_div tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_evaluates_on_edges tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_gradient tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_id tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_index tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_negation tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_printing tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_r_average tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_sign tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_unary_operator tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_x_average tests/unit/test_expression_tree/test_unary_operators.py::TestUnaryOperators::test_yz_average tests/unit/test_spatial_methods/test_base_spatial_method.py::TestSpatialMethod::test_broadcast_checks tests/unit/test_spatial_methods/test_base_spatial_method.py::TestSpatialMethod::test_discretise_spatial_variable tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_boundary_value_domain tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_concatenation tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_definite_integral tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_definite_integral_vector tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_delta_function tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_discretise_diffusivity_times_spatial_operator tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_discretise_spatial_variable tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_grad_1plus1d tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_and_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_grad_div_shapes_mixed_domain tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_grad_div_with_bcs_on_tab tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_heaviside tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_indefinite_integral tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_indefinite_integral_of_broadcasted_to_cell_edges tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_indefinite_integral_on_nodes tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_jacobian tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_mass_matrix_shape tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_neg_pos_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_node_to_edge_to_node tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_p2d_mass_matrix_shape tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_p2d_spherical_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_p2d_spherical_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Dirichlet_bcs tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py::TestFiniteVolume::test_spherical_grad_div_shapes_Neumann_bcs tests/unit/test_spatial_methods/test_zero_dimensional_method.py::TestZeroDimensionalSpatialMethod::test_broadcast_checks tests/unit/test_spatial_methods/test_zero_dimensional_method.py::TestZeroDimensionalSpatialMethod::test_discretise_spatial_variable
pybamm/expression_tree/unary_operators.py pybamm/spatial_methods/finite_volume.py pybamm/discretisations/discretisation.py pybamm/spatial_methods/zero_dimensional_method.py CHANGELOG.md pybamm/spatial_methods/spatial_method.py
pybamm/expression_tree/unary_operators.py pybamm/spatial_methods/finite_volume.py pybamm/discretisations/discretisation.py pybamm/spatial_methods/zero_dimensional_method.py CHANGELOG.md pybamm/spatial_methods/spatial_method.py