diff --git a/interfaces/cython/cantera/examples/surface_chemistry/lithium_ion_battery.py b/interfaces/cython/cantera/examples/surface_chemistry/lithium_ion_battery.py
new file mode 100644
index 0000000000..bcc1bd6a55
--- /dev/null
+++ b/interfaces/cython/cantera/examples/surface_chemistry/lithium_ion_battery.py
@@ -0,0 +1,172 @@
+"""
+This example calculates the cell voltage of a lithium-ion battery at
+given temperature, pressure, current, and range of state of charge (SOC).
+
+The thermodynamics are based on a graphite anode and a LiCoO2 cathode,
+modeled using the 'BinarySolutionTabulatedThermo' class.
+Further required cell parameters are the electrolyte ionic resistance, the
+stoichiometry ranges of the active materials (electrode balancing), and the
+surface area of the active materials.
+
+The functionality of this example is presented in greater detail in a jupyter
+notebook as well as the reference (which also describes the derivation of the
+'BinarySolutionTabulatedThermo' class):
+
+Reference:
+M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation
+of the thermodynamics of lithium-ion battery intercalation materials in the
+open-source software Cantera,” Electrochim. Acta 323, 134797 (2019),
+https://doi.org/10.1016/j.electacta.2019.134797
+
+Requires: cantera >= 2.6.0, matplotlib >= 2.0
+Keywords: surface chemistry, kinetics, electrochemistry, battery, plotting
+"""
+
+import cantera as ct
+import numpy as np
+
+# Parameters
+samples = 101
+soc = np.linspace(0., 1., samples)  # [-] Input state of charge (0...1)
+current = -1  # [A] Externally-applied current, negative for discharge
+T = 293  # [K] Temperature
+P = ct.one_atm  # [Pa] Pressure
+
+# Cell properties
+input_file = "lithium_ion_battery.yaml"  # Cantera input file name
+R_electrolyte = 0.0384  # [Ohm] Electrolyte resistance
+area_cathode = 1.1167  # [m^2] Cathode total active material surface area
+area_anode = 0.7824  # [m^2] Anode total active material surface area
+
+# Electrode balancing: The "balancing" of the electrodes relates the chemical
+# composition (lithium mole fraction in the active materials) to the macroscopic
+# cell-level state of charge.
+X_Li_anode_0 = 0.01  # [-] anode Li mole fraction at SOC = 0
+X_Li_anode_1 = 0.75  # [-] anode Li mole fraction at SOC = 100
+X_Li_cathode_0 = 0.99  # [-] cathode Li mole fraction at SOC = 0
+X_Li_cathode_1 = 0.49  # [-] cathode Li mole fraction at SOC = 100
+
+# Calculate mole fractions from SOC
+X_Li_anode = (X_Li_anode_1 - X_Li_anode_0) * soc + X_Li_anode_0
+X_Li_cathode = (X_Li_cathode_0 - X_Li_cathode_1) * (1 - soc) + X_Li_cathode_1
+
+# Import all Cantera phases
+anode = ct.Solution(input_file, "anode")
+cathode = ct.Solution(input_file, "cathode")
+metal = ct.Solution(input_file, "electron")
+electrolyte = ct.Solution(input_file, "electrolyte")
+anode_int = ct.Interface(
+    input_file, "edge_anode_electrolyte", adjacent=[anode, metal, electrolyte])
+cathode_int = ct.Interface(
+    input_file, "edge_cathode_electrolyte", adjacent=[cathode, metal, electrolyte])
+
+# Set the temperatures and pressures of all phases
+for phase in [anode, cathode, metal, electrolyte, anode_int, cathode_int]:
+    phase.TP = T, P
+
+
+# Root finding function
+def newton_solve(f, xstart, C=0.0):
+    """
+    Solve f(x) = C by Newton iteration using initial guess *xstart*
+    """
+    f0 = f(xstart) - C
+    x0 = xstart
+    dx = 1.0e-6
+    n = 0
+    while n < 200:
+        ff = f(x0 + dx) - C
+        dfdx = (ff - f0) / dx
+        step = - f0 / dfdx
+
+        # avoid taking steps too large
+        if abs(step) > 0.1:
+            step = 0.1 * step / abs(step)
+
+        x0 += step
+        emax = 0.00001  # 0.01 mV tolerance
+        if abs(f0) < emax and n > 8:
+            return x0
+        f0 = f(x0) - C
+        n += 1
+    raise Exception("no root!")
+
+
+# This function returns the Cantera calculated anode current (in A)
+def anode_current(phi_s, phi_l, X_Li_anode):
+    """
+    Current from the anode as a function of anode potential relative to
+    electrolyte.
+    """
+    # Set the active material mole fraction
+    anode.X = {"Li[anode]": X_Li_anode, "V[anode]": 1 - X_Li_anode}
+
+    # Set the electrode and electrolyte potential
+    metal.electric_potential = phi_s
+    electrolyte.electric_potential = phi_l
+
+    # Get the net reaction rate at the anode-side interface
+    # Reaction according to input file:
+    # Li+[electrolyte] + V[anode] + electron <=> Li[anode]
+    r = anode_int.net_rates_of_progress  # [kmol/m2/s]
+
+    # Calculate the current. Should be negative for cell discharge.
+    return r * ct.faraday * area_anode
+
+
+# This function returns the Cantera calculated cathode current (in A)
+def cathode_current(phi_s, phi_l, X_Li_cathode):
+    """
+    Current to the cathode as a function of cathode potential relative to electrolyte
+    """
+    # Set the active material mole fractions
+    cathode.X = {"Li[cathode]": X_Li_cathode, "V[cathode]": 1 - X_Li_cathode}
+
+    # Set the electrode and electrolyte potential
+    metal.electric_potential = phi_s
+    electrolyte.electric_potential = phi_l
+
+    # Get the net reaction rate at the cathode-side interface
+    # Reaction according to input file:
+    # Li+[electrolyte] + V[cathode] + electron <=> Li[cathode]
+    r = cathode_int.net_rates_of_progress  # [kmol/m2/s]
+
+    # Calculate the current. Should be negative for cell discharge.
+    return - r * ct.faraday * area_cathode
+
+
+# Calculate cell voltage, separately for each entry of the input vectors
+V_cell = np.zeros_like(soc)
+phi_l_anode = 0
+phi_s_cathode = 0
+for i in range(samples):
+    # Set anode electrode potential to 0
+    phi_s_anode = 0
+
+    # Calculate anode electrolyte potential
+    phi_l_anode = newton_solve(
+        lambda E: anode_current(phi_s_anode, E, X_Li_anode[i]),
+        phi_l_anode, C=current)
+
+    # Calculate cathode electrolyte potential
+    phi_l_cathode = phi_l_anode + current * R_electrolyte
+
+    # Calculate cathode electrode potential
+    phi_s_cathode = newton_solve(
+        lambda E: cathode_current(E, phi_l_cathode, X_Li_cathode[i]),
+        phi_s_cathode, C=current)
+
+    # Calculate cell voltage
+    V_cell[i] = phi_s_cathode - phi_s_anode
+
+try:
+    import matplotlib.pyplot as plt
+
+    # Plot the cell voltage, as a function of the state of charge
+    plt.plot(soc * 100, V_cell)
+    plt.xlabel("State of charge / %")
+    plt.ylabel("Cell voltage / V")
+    plt.show()
+
+except ImportError:
+    print("Install matplotlib to plot the outputs")
diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py
index 7c1cb7b5fd..c9e9113427 100644
--- a/interfaces/cython/cantera/test/test_kinetics.py
+++ b/interfaces/cython/cantera/test/test_kinetics.py
@@ -916,6 +916,106 @@ class TestSofcKinetics3(TestSofcKinetics):
     _mech = "sofc.cti"
 
 
+class TestLithiumIonBatteryKinetics(utilities.CanteraTest):
+    """ Test based on lithium_ion_battery.py """
+    _mech = "lithium_ion_battery.yaml"
+
+    def test_lithium_ion_battery(self):
+        mech = self._mech
+        samples = 11
+        soc = np.linspace(0., 1., samples)  # [-] Input state of charge (0...1)
+        current = -1  # [A] Externally-applied current, negative for discharge
+        T = 293  # T in K
+        P = ct.one_atm
+        R_electrolyte = 0.0384  # [Ohm] Electrolyte resistance
+        area_cathode = 1.1167  # [m^2] Cathode total active material surface area
+        area_anode = 0.7824  # [m^2] Anode total active material surface area
+
+        # Calculate mole fractions from SOC
+        X_Li_anode_0 = 0.01  # [-] anode Li mole fraction at SOC = 0
+        X_Li_anode_1 = 0.75  # [-] anode Li mole fraction at SOC = 100
+        X_Li_cathode_0 = 0.99  # [-] cathode Li mole fraction at SOC = 0
+        X_Li_cathode_1 = 0.49  # [-] cathode Li mole fraction at SOC = 100
+        X_Li_anode = (X_Li_anode_1 - X_Li_anode_0) * soc + X_Li_anode_0
+        X_Li_cathode = (X_Li_cathode_0 - X_Li_cathode_1) * (1 - soc) + X_Li_cathode_1
+
+        def newton_solve(f, xstart, C=0.0):
+            """ Solve f(x) = C by Newton iteration. """
+            x0 = xstart
+            dx = 1.0e-6
+            n = 0
+            while True:
+                n += 1
+                f0 = f(x0) - C
+                x0 -= f0 / (f(x0 + dx) - C - f0) * dx
+                if n > 1000:
+                    raise Exception('No convergence in Newton solve')
+                if abs(f0) < 0.00001:
+                    return x0
+
+        # Phases
+        anode, cathode, metal, electrolyte = ct.import_phases(
+            mech, ["anode", "cathode", "electron", "electrolyte"])
+        anode_int = ct.Interface(
+            mech, "edge_anode_electrolyte", adjacent=[anode, metal, electrolyte])
+        cathode_int = ct.Interface(
+            mech, "edge_cathode_electrolyte", adjacent=[cathode, metal, electrolyte])
+
+        # initialization
+        for phase in [anode, cathode, metal, electrolyte, anode_int, cathode_int]:
+            phase.TP = T, P
+
+        # This function returns the Cantera calculated anode current
+        def anode_current(phi_s, phi_l, X_Li_anode):
+            # Set mole fraction and electrode and electrolyte potential
+            anode.X = {"Li[anode]": X_Li_anode, "V[anode]": 1 - X_Li_anode}
+            metal.electric_potential = phi_s
+            electrolyte.electric_potential = phi_l
+
+            # Calculate the current.
+            return ct.faraday * anode_int.net_rates_of_progress * area_anode
+
+        # This function returns the Cantera calculated cathode current
+        def cathode_current(phi_s, phi_l, X_Li_cathode):
+            # Set mole fraction and electrode and electrolyte potential
+            cathode.X = {"Li[cathode]": X_Li_cathode, "V[cathode]": 1 - X_Li_cathode}
+            metal.electric_potential = phi_s
+            electrolyte.electric_potential = phi_l
+
+            # Calculate the current. Should be negative for cell discharge.
+            return - ct.faraday * cathode_int.net_rates_of_progress * area_cathode
+
+        # Calculate cell voltage, separately for each entry of the input vectors
+        data = []
+        phi_l_anode = 0
+        phi_s_cathode = 0
+        for i in range(samples):
+            # Calculate anode electrolyte potential
+            phi_s_anode = 0
+            phi_l_anode = newton_solve(
+                lambda E: anode_current(phi_s_anode, E, X_Li_anode[i]),
+                phi_l_anode, C=current)
+
+            # Calculate cathode electrode potential
+            phi_l_cathode = phi_l_anode + current * R_electrolyte
+            phi_s_cathode = newton_solve(
+                lambda E: cathode_current(E, phi_l_cathode, X_Li_cathode[i]),
+                phi_s_cathode, C=current)
+
+            # Calculate cell voltage
+            data.append(phi_s_cathode - phi_s_anode)
+
+        data = np.array(data).ravel()
+        ref = np.genfromtxt(self.test_data_path / "lithium-ion-battery-test.csv")
+        assert np.allclose(data, ref, rtol=1e-7)
+
+
+@pytest.mark.usefixtures("allow_deprecated")
+class TestLithiumIonBatteryKinetics2(TestLithiumIonBatteryKinetics):
+    """ Test using legacy framework; included to retain coverage """
+    _mech = "lithium_ion_battery.cti"
+
+
 class TestDuplicateReactions(utilities.CanteraTest):
     infile = 'duplicate-reactions.yaml'
 
diff --git a/test/data/lithium-ion-battery-test.csv b/test/data/lithium-ion-battery-test.csv
new file mode 100644
index 0000000000..b9ed0b778c
--- /dev/null
+++ b/test/data/lithium-ion-battery-test.csv
@@ -0,0 +1,11 @@
+2.65335507
+3.57418843
+3.64605784
+3.68174416
+3.70084855
+3.71586159
+3.75736585
+3.81434722
+3.868883
+3.91854252
+4.01938813