diff --git a/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst b/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst
index 1c7b8b5049..e1876dfc57 100644
--- a/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst
+++ b/docs/source/models/submodels/external_circuit/function_control_external_circuit.rst
@@ -8,4 +8,7 @@ Function control external circuit
     :members:
 
 .. autoclass:: pybamm.external_circuit.PowerFunctionControl
+    :members:
+
+.. autoclass:: pybamm.external_circuit.CCCVFunctionControl
     :members:
\ No newline at end of file
diff --git a/pybamm/experiments/experiment.py b/pybamm/experiments/experiment.py
index 676690231e..f0fcbb60d3 100644
--- a/pybamm/experiments/experiment.py
+++ b/pybamm/experiments/experiment.py
@@ -55,7 +55,11 @@ class Experiment:
         faster at simulating individual steps but has higher set-up overhead
     drive_cycles : dict
         Dictionary of drive cycles to use for this experiment.
-
+    cccv_handling : str, optional
+        How to handle CCCV. If "two-step" (default), then the experiment is run in
+        two steps (CC then CV). If "ode", then the experiment is run in a single step
+        using an ODE for current: see
+        :class:`pybamm.external_circuit.CCCVFunctionControl` for details.
     """
 
     def __init__(
@@ -66,7 +70,11 @@ def __init__(
         termination=None,
         use_simulation_setup_type="new",
         drive_cycles={},
+        cccv_handling="two-step",
     ):
+        if cccv_handling not in ["two-step", "ode"]:
+            raise ValueError("cccv_handling should be either 'two-step' or 'ode'")
+        self.cccv_handling = cccv_handling
 
         self.period = self.convert_time_to_seconds(period.split())
         operating_conditions_cycles = []
@@ -75,9 +83,28 @@ def __init__(
             if (isinstance(cycle, tuple) or isinstance(cycle, str)) and all(
                 [isinstance(cond, str) for cond in cycle]
             ):
-                operating_conditions_cycles.append(
-                    cycle if isinstance(cycle, tuple) else (cycle,)
-                )
+                if isinstance(cycle, str):
+                    processed_cycle = (cycle,)
+                else:
+                    processed_cycle = []
+                    idx = 0
+                    finished = False
+                    while not finished:
+                        step = cycle[idx]
+                        if idx < len(cycle) - 1:
+                            next_step = cycle[idx + 1]
+                        else:
+                            next_step = None
+                            finished = True
+                        if self.is_cccv(step, next_step):
+                            processed_cycle.append(step + " then " + next_step)
+                            idx += 2
+                        else:
+                            processed_cycle.append(step)
+                            idx += 1
+                        if idx >= len(cycle):
+                            finished = True
+                operating_conditions_cycles.append(tuple(processed_cycle))
             else:
                 try:
                     # Condition is not a string
@@ -153,7 +180,20 @@ def read_string(self, cond, drive_cycles):
             must be numbers, 'C' denotes the unit of the external circuit (can be A for
             current, C for C-rate, V for voltage or W for power), and 'hours' denotes
             the unit of time (can be second(s), minute(s) or hour(s))
+        drive_cycles: dict
+            A map specifying the drive cycles
         """
+        if " then " in cond:
+            # If the string contains " then ", then this is a two-step CCCV experiment
+            # and we need to split it into two strings
+            cond_CC, cond_CV = cond.split(" then ")
+            op_CC, _ = self.read_string(cond_CC, drive_cycles)
+            op_CV, event_CV = self.read_string(cond_CV, drive_cycles)
+            return {
+                "electric": op_CC["electric"] + op_CV["electric"],
+                "time": op_CV["time"],
+                "period": op_CV["period"],
+            }, event_CV
         # Read period
         if " period)" in cond:
             cond, time_period = cond.split("(")
@@ -165,20 +205,12 @@ def read_string(self, cond, drive_cycles):
         if "Run" in cond:
             cond_list = cond.split()
             if "at" in cond:
-                raise ValueError(
-                    """Instruction must be
-                    For example: {}""".format(
-                        examples
-                    )
-                )
+                raise ValueError(f"Instruction must be of the form: {examples}")
             dc_types = ["(A)", "(V)", "(W)"]
             if all(x not in cond for x in dc_types):
                 raise ValueError(
-                    """Type of drive cycle must be
-                    specified using '(A)', '(V)' or '(W)'.
-                    For example: {}""".format(
-                        examples
-                    )
+                    "Type of drive cycle must be specified using '(A)', '(V)' or '(W)'."
+                    f" For example: {examples}"
                 )
             # Check for Events
             elif "for" in cond:
@@ -208,7 +240,7 @@ def read_string(self, cond, drive_cycles):
                 time = drive_cycles[cond_list[1]][:, 0][-1]
                 period = np.min(np.diff(drive_cycles[cond_list[1]][:, 0]))
                 events = None
-        elif "Run" not in cond:
+        else:
             if "for" in cond and "or until" in cond:
                 # e.g. for 3 hours or until 4.2 V
                 cond_list = cond.split()
@@ -238,7 +270,8 @@ def read_string(self, cond, drive_cycles):
                         examples
                     )
                 )
-        return electric + (time,) + (period,), events
+
+        return {"electric": electric, "time": time, "period": period}, events
 
     def extend_drive_cycle(self, drive_cycle, end_time):
         "Extends the drive cycle to enable for event"
@@ -404,3 +437,25 @@ def read_termination(self, termination):
                     "e.g. '80% capacity' or '4 Ah capacity'"
                 )
         return termination_dict
+
+    def is_cccv(self, step, next_step):
+        """
+        Check whether a step and the next step indicate a CCCV charge
+        """
+        if self.cccv_handling == "two-step" or next_step is None:
+            return False
+        # e.g. step="Charge at 2.0 A until 4.2V"
+        # next_step="Hold at 4.2V until C/50"
+        if (
+            step.startswith("Charge")
+            and "until" in step
+            and "V" in step
+            and "Hold at " in next_step
+            and "V until" in next_step
+        ):
+            _, events = self.read_string(step, None)
+            next_op, _ = self.read_string(next_step, None)
+            # Check that the event conditions are the same as the hold conditions
+            if events == next_op["electric"]:
+                return True
+        return False
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index a6e23efeb4..db6fc3a543 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -164,7 +164,7 @@ def __init__(self, extra_options):
             "lithium plating": ["none", "reversible", "irreversible"],
             "lithium plating porosity change": ["true", "false"],
             "loss of active material": ["none", "stress-driven", "reaction-driven"],
-            "operating mode": ["current", "voltage", "power"],
+            "operating mode": ["current", "voltage", "power", "CCCV"],
             "particle": [
                 "Fickian diffusion",
                 "fast diffusion",
@@ -410,7 +410,7 @@ def default_parameter_values(self):
     def default_geometry(self):
         return pybamm.battery_geometry(
             options=self.options,
-            current_collector_dimension=self.options["dimensionality"]
+            current_collector_dimension=self.options["dimensionality"],
         )
 
     @property
@@ -440,12 +440,8 @@ def default_submesh_types(self):
             "positive electrode": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
             "negative particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
             "positive particle": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
-            "negative particle size": pybamm.MeshGenerator(
-                pybamm.Uniform1DSubMesh
-            ),
-            "positive particle size": pybamm.MeshGenerator(
-                pybamm.Uniform1DSubMesh
-            ),
+            "negative particle size": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
+            "positive particle size": pybamm.MeshGenerator(pybamm.Uniform1DSubMesh),
         }
         if self.options["dimensionality"] == 0:
             base_submeshes["current collector"] = pybamm.MeshGenerator(pybamm.SubMesh0D)
@@ -731,7 +727,7 @@ def build_model(self):
         pybamm.logger.info("Finish building {}".format(self.name))
 
     def new_empty_copy(self):
-        """ See :meth:`pybamm.BaseModel.new_empty_copy()` """
+        """See :meth:`pybamm.BaseModel.new_empty_copy()`"""
         new_model = self.__class__(name=self.name, options=self.options)
         new_model.use_jacobian = self.use_jacobian
         new_model.convert_to_format = self.convert_to_format
@@ -756,6 +752,10 @@ def set_external_circuit_submodel(self):
             self.submodels[
                 "external circuit"
             ] = pybamm.external_circuit.PowerFunctionControl(self.param)
+        elif self.options["operating mode"] == "CCCV":
+            self.submodels[
+                "external circuit"
+            ] = pybamm.external_circuit.CCCVFunctionControl(self.param)
         elif callable(self.options["operating mode"]):
             self.submodels[
                 "external circuit"
diff --git a/pybamm/models/submodels/external_circuit/__init__.py b/pybamm/models/submodels/external_circuit/__init__.py
index f181a9cca9..6cc5b14220 100644
--- a/pybamm/models/submodels/external_circuit/__init__.py
+++ b/pybamm/models/submodels/external_circuit/__init__.py
@@ -4,6 +4,7 @@
     FunctionControl,
     VoltageFunctionControl,
     PowerFunctionControl,
+    CCCVFunctionControl,
     LeadingOrderFunctionControl,
     LeadingOrderVoltageFunctionControl,
     LeadingOrderPowerFunctionControl,
diff --git a/pybamm/models/submodels/external_circuit/current_control_external_circuit.py b/pybamm/models/submodels/external_circuit/current_control_external_circuit.py
index 111efa3a19..c51bab2fbb 100644
--- a/pybamm/models/submodels/external_circuit/current_control_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/current_control_external_circuit.py
@@ -1,6 +1,7 @@
 #
 # External circuit with current control
 #
+import pybamm
 from .base_external_circuit import BaseModel, LeadingOrderBaseModel
 
 
@@ -17,6 +18,7 @@ def get_fundamental_variables(self):
         I = self.param.dimensional_current_with_time
 
         variables = {
+            "Current density variable": pybamm.Scalar(1, name="i_cell"),
             "Total current density": i_cell,
             "Total current density [A.m-2]": i_cell_dim,
             "Current [A]": I,
diff --git a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
index aa53e23c86..ad1878beb6 100644
--- a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
@@ -6,22 +6,41 @@
 
 
 class FunctionControl(BaseModel):
-    """External circuit with an arbitrary function."""
+    """
+    External circuit with an arbitrary function, implemented as a control on the current
+    either via an algebraic equation, or a differential equation.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    external_circuit_function : callable
+        The function that controls the current
+    control : str, optional
+        The type of control to use. Must be one of 'algebraic' (default)
+        or 'differential'.
+    """
 
-    def __init__(self, param, external_circuit_function):
+    def __init__(self, param, external_circuit_function, control="algebraic"):
         super().__init__(param)
         self.external_circuit_function = external_circuit_function
+        self.control = control
 
     def get_fundamental_variables(self):
         param = self.param
         # Current is a variable
-        i_cell = pybamm.Variable("Total current density")
+        i_var = pybamm.Variable("Current density variable")
+        if self.control == "algebraic":
+            i_cell = i_var
+        elif self.control == "differential":
+            i_cell = pybamm.maximum(i_var, param.current_with_time)
 
         # Update derived variables
-        I = i_cell * pybamm.AbsoluteValue(param.I_typ)
+        I = i_cell * abs(param.I_typ)
         i_cell_dim = I / (param.n_electrodes_parallel * param.A_cc)
 
         variables = {
+            "Current density variable": i_var,
             "Total current density": i_cell,
             "Total current density [A.m-2]": i_cell_dim,
             "Current [A]": I,
@@ -36,15 +55,25 @@ def get_fundamental_variables(self):
     def set_initial_conditions(self, variables):
         super().set_initial_conditions(variables)
         # Initial condition as a guess for consistent initial conditions
-        i_cell = variables["Total current density"]
+        i_cell = variables["Current density variable"]
         self.initial_conditions[i_cell] = self.param.current_with_time
 
+    def set_rhs(self, variables):
+        super().set_rhs(variables)
+        # External circuit submodels are always equations on the current
+        # The external circuit function should provide an update law for the current
+        # based on current/voltage/power/etc.
+        if self.control == "differential":
+            i_cell = variables["Current density variable"]
+            self.rhs[i_cell] = self.external_circuit_function(variables)
+
     def set_algebraic(self, variables):
         # External circuit submodels are always equations on the current
         # The external circuit function should fix either the current, or the voltage,
         # or a combination (e.g. I*V for power control)
-        i_cell = variables["Total current density"]
-        self.algebraic[i_cell] = self.external_circuit_function(variables)
+        if self.control == "algebraic":
+            i_cell = variables["Current density variable"]
+            self.algebraic[i_cell] = self.external_circuit_function(variables)
 
 
 class VoltageFunctionControl(FunctionControl):
@@ -53,7 +82,7 @@ class VoltageFunctionControl(FunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, self.constant_voltage)
+        super().__init__(param, self.constant_voltage, control="algebraic")
 
     def constant_voltage(self, variables):
         V = variables["Terminal voltage [V]"]
@@ -66,7 +95,7 @@ class PowerFunctionControl(FunctionControl):
     """External circuit with power control."""
 
     def __init__(self, param):
-        super().__init__(param, self.constant_power)
+        super().__init__(param, self.constant_power, control="algebraic")
 
     def constant_power(self, variables):
         I = variables["Current [A]"]
@@ -76,11 +105,29 @@ def constant_power(self, variables):
         )
 
 
+class CCCVFunctionControl(FunctionControl):
+    """External circuit with constant-current constant-voltage control."""
+
+    def __init__(self, param):
+        super().__init__(param, self.cccv, control="differential")
+
+    def cccv(self, variables):
+        # Multiply by the time scale so that the votage overshoot only lasts a few
+        # seconds
+        K_aw = 1 * self.param.timescale  # anti-windup
+        K_V = 1 * self.param.timescale
+        i_var = variables["Current density variable"]
+        i_cell = variables["Total current density"]
+        V = variables["Terminal voltage [V]"]
+        V_CCCV = pybamm.Parameter("Voltage function [V]")
+        return -K_aw * (i_var - i_cell) + K_V * (V - V_CCCV)
+
+
 class LeadingOrderFunctionControl(FunctionControl, LeadingOrderBaseModel):
     """External circuit with an arbitrary function, at leading order."""
 
-    def __init__(self, param, external_circuit_class):
-        super().__init__(param, external_circuit_class)
+    def __init__(self, param, external_circuit_class, control="algebraic"):
+        super().__init__(param, external_circuit_class, control=control)
 
     def _get_current_variable(self):
         return pybamm.Variable("Leading-order total current density")
@@ -93,7 +140,7 @@ class LeadingOrderVoltageFunctionControl(LeadingOrderFunctionControl):
     """
 
     def __init__(self, param):
-        super().__init__(param, self.constant_voltage)
+        super().__init__(param, self.constant_voltage, control="algebraic")
 
     def constant_voltage(self, variables):
         V = variables["Terminal voltage [V]"]
@@ -106,7 +153,7 @@ class LeadingOrderPowerFunctionControl(LeadingOrderFunctionControl):
     """External circuit with power control, at leading order."""
 
     def __init__(self, param):
-        super().__init__(param, self.constant_power)
+        super().__init__(param, self.constant_power, control="algebraic")
 
     def constant_power(self, variables):
         I = variables["Current [A]"]
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 52baca3646..000960e3dc 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -41,19 +41,6 @@ def constant_current_constant_voltage_constant_power(variables):
     )
 
 
-def constant_voltage(variables, V_applied):
-    V = variables["Terminal voltage [V]"]
-    n_cells = pybamm.Parameter("Number of cells connected in series to make a battery")
-    return V - V_applied / n_cells
-
-
-def constant_power(variables, P_applied):
-    I = variables["Current [A]"]
-    V = variables["Terminal voltage [V]"]
-    n_cells = pybamm.Parameter("Number of cells connected in series to make a battery")
-    return V * I - P_applied / n_cells
-
-
 class Simulation:
     """A Simulation class for easy building and running of PyBaMM simulations.
 
@@ -101,7 +88,8 @@ def __init__(
             if experiment is not None:
                 raise NotImplementedError(
                     "BasicDFNHalfCell is not compatible "
-                    "with experiment simulations yet.")
+                    "with experiment simulations yet."
+                )
 
         if experiment is None:
             # Check to see if the current is provided as data (i.e. drive cycle)
@@ -158,6 +146,8 @@ def set_up_experiment(self, model, experiment):
         operating condition in the experiment. The model will then be solved by
         integrating the model successively with each group of inputs, one group at a
         time.
+        This needs to be done here and not in the Experiment class because the nominal
+        cell capacity (from the parameters) is used to convert C-rate to current.
         """
         self.operating_mode = "with experiment"
 
@@ -173,48 +163,51 @@ def set_up_experiment(self, model, experiment):
         self._experiment_inputs = []
         self._experiment_times = []
         for op, events in zip(experiment.operating_conditions, experiment.events):
-            if op[1] in ["A", "C"]:
-                # Update inputs for constant current
+            operating_inputs = {
+                "Current switch": 0,
+                "Voltage switch": 0,
+                "Power switch": 0,
+                "CCCV switch": 0,
+                "Current input [A]": 0,
+                "Voltage input [V]": 0,  # doesn't matter
+                "Power input [W]": 0,  # doesn't matter
+            }
+            op_control = op["electric"][1]
+            if op_control in ["A", "C"]:
                 capacity = self._parameter_values["Nominal cell capacity [A.h]"]
-                if op[1] == "A":
-                    I = op[0]
+                if op_control == "A":
+                    I = op["electric"][0]
                     Crate = I / capacity
                 else:
                     # Scale C-rate with capacity to obtain current
-                    Crate = op[0]
+                    Crate = op["electric"][0]
                     I = Crate * capacity
-                operating_inputs = {
-                    "Current switch": 1,
-                    "Voltage switch": 0,
-                    "Power switch": 0,
-                    "Current input [A]": I,
-                    "Voltage input [V]": 0,  # doesn't matter
-                    "Power input [W]": 0,  # doesn't matter
-                }
-            elif op[1] == "V":
+                if len(op["electric"]) == 4:
+                    # Update inputs for CCCV
+                    op_control = "CCCV"  # change to CCCV
+                    V = op["electric"][2]
+                    operating_inputs.update(
+                        {
+                            "CCCV switch": 1,
+                            "Current input [A]": I,
+                            "Voltage input [V]": V,
+                        }
+                    )
+                else:
+                    # Update inputs for constant current
+                    operating_inputs.update(
+                        {"Current switch": 1, "Current input [A]": I}
+                    )
+            elif op_control == "V":
                 # Update inputs for constant voltage
-                V = op[0]
-                operating_inputs = {
-                    "Current switch": 0,
-                    "Voltage switch": 1,
-                    "Power switch": 0,
-                    "Current input [A]": 0,  # doesn't matter
-                    "Voltage input [V]": V,
-                    "Power input [W]": 0,  # doesn't matter
-                }
-            elif op[1] == "W":
+                V = op["electric"][0]
+                operating_inputs.update({"Voltage switch": 1, "Voltage input [V]": V})
+            elif op_control == "W":
                 # Update inputs for constant power
-                P = op[0]
-                operating_inputs = {
-                    "Current switch": 0,
-                    "Voltage switch": 0,
-                    "Power switch": 1,
-                    "Current input [A]": 0,  # doesn't matter
-                    "Voltage input [V]": 0,  # doesn't matter
-                    "Power input [W]": P,
-                }
+                P = op["electric"][0]
+                operating_inputs.update({"Power switch": 1, "Power input [W]": P})
             # Update period
-            operating_inputs["period"] = op[3]
+            operating_inputs["period"] = op["period"]
             # Update events
             if events is None:
                 # make current and voltage values that won't be hit
@@ -241,12 +234,14 @@ def set_up_experiment(self, model, experiment):
 
             self._experiment_inputs.append(operating_inputs)
             # Add time to the experiment times
-            dt = op[2]
+            dt = op["time"]
             if dt is None:
-                if op[1] in ["A", "C"]:
+                if op_control in ["A", "C", "CCCV"]:
                     # Current control: max simulation time: 3 * max simulation time
                     # based on C-rate
                     dt = 3 / abs(Crate) * 3600  # seconds
+                    if op_control == "CCCV":
+                        dt *= 5  # 5x longer for CCCV
                 else:
                     # max simulation time: 1 day
                     dt = 24 * 3600  # seconds
@@ -328,9 +323,13 @@ def set_up_model_for_experiment_old(self, model):
 
         self.model = new_model
 
+        operating_conditions = set(
+            x["electric"] + (x["time"],) + (x["period"],)
+            for x in self.experiment.operating_conditions
+        )
         self.op_conds_to_model_and_param = {
             op_cond[:2]: (new_model, self.parameter_values)
-            for op_cond in set(self.experiment.operating_conditions)
+            for op_cond in operating_conditions
         }
 
     def set_up_model_for_experiment_new(self, model):
@@ -347,7 +346,7 @@ def set_up_model_for_experiment_new(self, model):
         ):
             # Create model for this operating condition if it has not already been seen
             # before
-            if op_cond[:2] not in self.op_conds_to_model_and_param:
+            if op_cond["electric"] not in self.op_conds_to_model_and_param:
                 if op_inputs["Current switch"] == 1:
                     # Current control
                     # Make a new copy of the model (we will update events later))
@@ -358,10 +357,16 @@ def set_up_model_for_experiment_new(self, model):
                     # To do so, we replace all instances of the current density in the
                     # model with a current density variable, which is obtained from the
                     # FunctionControl submodel
+                    # check which kind of external circuit model we need (differential
+                    # or algebraic)
+                    if op_inputs["CCCV switch"] == 1:
+                        control = "differential"
+                    else:
+                        control = "algebraic"
                     # create the FunctionControl submodel and extract variables
                     external_circuit_variables = (
                         pybamm.external_circuit.FunctionControl(
-                            model.param, None
+                            model.param, None, control=control
                         ).get_fundamental_variables()
                     )
 
@@ -373,15 +378,15 @@ def set_up_model_for_experiment_new(self, model):
                     replacer = pybamm.SymbolReplacer(symbol_replacement_map)
                     new_model = replacer.process_model(model, inplace=False)
 
-                    # Update the algebraic equation and initial conditions for
+                    # Update the rhs or algebraic equation and initial conditions for
                     # FunctionControl
-                    # This creates an algebraic equation for the current to allow
-                    # current, voltage, or power control, together with the appropriate
-                    # guess for the initial condition.
+                    # This creates a differential or algebraic equation for the current
+                    # to allow current, voltage, or power control, together with the
+                    # appropriate guess for the initial condition.
                     # External circuit submodels are always equations on the current
                     # The external circuit function should fix either the current, or
                     # the voltage, or a combination (e.g. I*V for power control)
-                    i_cell = new_model.variables["Total current density"]
+                    i_cell = new_model.variables["Current density variable"]
                     new_model.initial_conditions[
                         i_cell
                     ] = new_model.param.current_with_time
@@ -403,14 +408,28 @@ def set_up_model_for_experiment_new(self, model):
                         ]
                     )
                     if op_inputs["Voltage switch"] == 1:
-                        new_model.algebraic[i_cell] = constant_voltage(
-                            new_model.variables,
-                            pybamm.Parameter("Voltage function [V]"),
+                        new_model.algebraic[
+                            i_cell
+                        ] = pybamm.external_circuit.VoltageFunctionControl(
+                            new_model.param
+                        ).constant_voltage(
+                            new_model.variables
                         )
                     elif op_inputs["Power switch"] == 1:
-                        new_model.algebraic[i_cell] = constant_power(
-                            new_model.variables,
-                            pybamm.Parameter("Power function [W]"),
+                        new_model.algebraic[
+                            i_cell
+                        ] = pybamm.external_circuit.PowerFunctionControl(
+                            new_model.param
+                        ).constant_power(
+                            new_model.variables
+                        )
+                    elif op_inputs["CCCV switch"] == 1:
+                        new_model.algebraic[
+                            i_cell
+                        ] = pybamm.external_circuit.CCCVFunctionControl(
+                            new_model.param
+                        ).cccv(
+                            new_model.variables
                         )
 
                 # add voltage events to the model
@@ -449,8 +468,16 @@ def set_up_model_for_experiment_new(self, model):
                         {"Power function [W]": op_inputs["Power input [W]"]},
                         check_already_exists=False,
                     )
+                elif op_inputs["CCCV switch"] == 1:
+                    new_parameter_values.update(
+                        {
+                            "Current function [A]": op_inputs["Current input [A]"],
+                            "Voltage function [V]": op_inputs["Voltage input [V]"],
+                        },
+                        check_already_exists=False,
+                    )
 
-                self.op_conds_to_model_and_param[op_cond[:2]] = (
+                self.op_conds_to_model_and_param[op_cond["electric"]] = (
                     new_model,
                     new_parameter_values,
                 )
@@ -627,7 +654,7 @@ def solve(
                 }
             )
             # For experiments also update the following
-            if hasattr(self, 'op_conds_to_model_and_param'):
+            if hasattr(self, "op_conds_to_model_and_param"):
                 for key, (model, param) in self.op_conds_to_model_and_param.items():
                     param.update(
                         {
@@ -782,7 +809,9 @@ def solve(
                     exp_inputs = self._experiment_inputs[idx]
                     dt = self._experiment_times[idx]
                     op_conds_str = self.experiment.operating_conditions_strings[idx]
-                    op_conds_elec = self.experiment.operating_conditions[idx][:2]
+                    op_conds_elec = self.experiment.operating_conditions[idx][
+                        "electric"
+                    ]
                     model = self.op_conds_to_built_models[op_conds_elec]
                     # Use 1-indexing for printing cycle number as it is more
                     # human-intuitive
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 46b68ed960..5971a45de1 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -990,7 +990,7 @@ def step(
                 model.y0 = old_solution.all_ys[-1][:, -1]
             else:
                 model.y0 = (
-                    model.set_initial_conditions_from(old_solution, inplace=False)
+                    model.set_initial_conditions_from(old_solution)
                     .concatenated_initial_conditions.evaluate(0, inputs=ext_and_inputs)
                     .flatten()
                 )
diff --git a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
index c3a2bc243a..fa6cb7ed05 100644
--- a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
+++ b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py
@@ -139,6 +139,29 @@ def constant_power(variables):
         I1 = solutions[1]["Current [A]"].entries
         np.testing.assert_array_equal(I0, I1)
 
+    def test_cccv(self):
+        # load models
+        model = pybamm.lithium_ion.SPM({"operating mode": "CCCV"})
+
+        # load parameter values and process models and geometry
+        param = model.default_parameter_values
+
+        # First model: 4W charge
+        param.update({"Voltage function [V]": 4.2}, check_already_exists=False)
+
+        # set parameters and discretise models
+        # create geometry
+        geometry = model.default_geometry
+        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)
+
+        # solve model
+        t_eval = np.linspace(0, 3600, 100)
+        model.default_solver.solve(model, t_eval)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_experiments/test_experiment.py b/tests/unit/test_experiments/test_experiment.py
index df88b4a937..db6f6196e2 100644
--- a/tests/unit/test_experiments/test_experiment.py
+++ b/tests/unit/test_experiments/test_experiment.py
@@ -2,7 +2,7 @@
 # Test the base experiment class
 #
 import pybamm
-import numpy
+import numpy as np
 import unittest
 import pandas as pd
 import os
@@ -43,51 +43,52 @@ def test_read_strings(self):
             period="20 seconds",
         )
 
+        self.assertEqual(
+            experiment.operating_conditions[:-3],
+            [
+                {"electric": (1, "C"), "time": 1800.0, "period": 20.0},
+                {"electric": (0.05, "C"), "time": 1800.0, "period": 20.0},
+                {"electric": (-0.5, "C"), "time": 2700.0, "period": 20.0},
+                {"electric": (1, "A"), "time": 1800.0, "period": 20.0},
+                {"electric": (-0.2, "A"), "time": 2700.0, "period": 60.0},
+                {"electric": (1, "W"), "time": 1800.0, "period": 20.0},
+                {"electric": (-0.2, "W"), "time": 2700.0, "period": 20.0},
+                {"electric": (0, "A"), "time": 600.0, "period": 300.0},
+                {"electric": (1, "V"), "time": 20.0, "period": 20.0},
+                {"electric": (-1, "C"), "time": None, "period": 20.0},
+                {"electric": (4.1, "V"), "time": None, "period": 20.0},
+                {"electric": (3, "V"), "time": None, "period": 20.0},
+                {"electric": (1 / 3, "C"), "time": 7200.0, "period": 20.0},
+            ],
+        )
         # Calculation for operating conditions of drive cycle
         time_0 = drive_cycle[:, 0][-1]
-        period_0 = numpy.min(numpy.diff(drive_cycle[:, 0]))
+        period_0 = np.min(np.diff(drive_cycle[:, 0]))
         drive_cycle_1 = experiment.extend_drive_cycle(drive_cycle, end_time=300)
         time_1 = drive_cycle_1[:, 0][-1]
-        period_1 = numpy.min(numpy.diff(drive_cycle_1[:, 0]))
+        period_1 = np.min(np.diff(drive_cycle_1[:, 0]))
         drive_cycle_2 = experiment.extend_drive_cycle(drive_cycle, end_time=1800)
         time_2 = drive_cycle_2[:, 0][-1]
-        period_2 = numpy.min(numpy.diff(drive_cycle_2[:, 0]))
-
-        self.assertEqual(
-            experiment.operating_conditions[:-3],
-            [
-                (1, "C", 1800.0, 20.0),
-                (0.05, "C", 1800.0, 20.0),
-                (-0.5, "C", 2700.0, 20.0),
-                (1, "A", 1800.0, 20.0),
-                (-0.2, "A", 2700.0, 60.0),
-                (1, "W", 1800.0, 20.0),
-                (-0.2, "W", 2700.0, 20.0),
-                (0, "A", 600.0, 300.0),
-                (1, "V", 20.0, 20.0),
-                (-1, "C", None, 20.0),
-                (4.1, "V", None, 20.0),
-                (3, "V", None, 20.0),
-                (1 / 3, "C", 7200.0, 20.0),
-            ],
-        )
+        period_2 = np.min(np.diff(drive_cycle_2[:, 0]))
         # Check drive cycle operating conditions
-        self.assertTrue(
-            (
-                (experiment.operating_conditions[-3][0] == drive_cycle).all()
-                & (experiment.operating_conditions[-3][1] == "A")
-                & (experiment.operating_conditions[-3][2] == time_0).all()
-                & (experiment.operating_conditions[-3][3] == period_0).all()
-                & (experiment.operating_conditions[-2][0] == drive_cycle_1).all()
-                & (experiment.operating_conditions[-2][1] == "V")
-                & (experiment.operating_conditions[-2][2] == time_1).all()
-                & (experiment.operating_conditions[-2][3] == period_1).all()
-                & (experiment.operating_conditions[-1][0] == drive_cycle_2).all()
-                & (experiment.operating_conditions[-1][1] == "W")
-                & (experiment.operating_conditions[-1][2] == time_2).all()
-                & (experiment.operating_conditions[-1][3] == period_2).all()
-            )
+        np.testing.assert_array_equal(
+            experiment.operating_conditions[-3]["electric"][0], drive_cycle
+        )
+        self.assertEqual(experiment.operating_conditions[-3]["electric"][1], "A")
+        self.assertEqual(experiment.operating_conditions[-3]["time"], time_0)
+        self.assertEqual(experiment.operating_conditions[-3]["period"], period_0)
+        np.testing.assert_array_equal(
+            experiment.operating_conditions[-2]["electric"][0], drive_cycle_1
+        )
+        self.assertEqual(experiment.operating_conditions[-2]["electric"][1], "V")
+        self.assertEqual(experiment.operating_conditions[-2]["time"], time_1)
+        self.assertEqual(experiment.operating_conditions[-2]["period"], period_1)
+        np.testing.assert_array_equal(
+            experiment.operating_conditions[-1]["electric"][0], drive_cycle_2
         )
+        self.assertEqual(experiment.operating_conditions[-1]["electric"][1], "W")
+        self.assertEqual(experiment.operating_conditions[-1]["time"], time_2)
+        self.assertEqual(experiment.operating_conditions[-1]["period"], period_2)
         self.assertEqual(
             experiment.events,
             [
@@ -112,6 +113,58 @@ def test_read_strings(self):
         self.assertEqual(experiment.parameters, {"test": "test"})
         self.assertEqual(experiment.period, 20)
 
+    def test_read_strings_cccv_combined(self):
+        experiment = pybamm.Experiment(
+            [
+                (
+                    "Discharge at C/20 for 0.5 hours",
+                    "Charge at 0.5 C until 1V",
+                    "Hold at 1V until C/50",
+                    "Discharge at C/20 for 0.5 hours",
+                ),
+            ],
+            cccv_handling="ode",
+        )
+        self.assertEqual(
+            experiment.operating_conditions,
+            [
+                {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0},
+                {"electric": (-0.5, "C", 1, "V"), "time": None, "period": 60.0},
+                {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0},
+            ],
+        )
+        self.assertEqual(experiment.events, [None, (0.02, "C"), None])
+
+        # Cases that don't quite match shouldn't do CCCV setup
+        experiment = pybamm.Experiment(
+            [
+                "Charge at 0.5 C until 2V",
+                "Hold at 1V until C/50",
+            ],
+            cccv_handling="ode",
+        )
+        self.assertEqual(
+            experiment.operating_conditions,
+            [
+                {"electric": (-0.5, "C"), "time": None, "period": 60.0},
+                {"electric": (1, "V"), "time": None, "period": 60.0},
+            ],
+        )
+        experiment = pybamm.Experiment(
+            [
+                "Charge at 0.5 C for 2 minutes",
+                "Hold at 1V until C/50",
+            ],
+            cccv_handling="ode",
+        )
+        self.assertEqual(
+            experiment.operating_conditions,
+            [
+                {"electric": (-0.5, "C"), "time": 120.0, "period": 60.0},
+                {"electric": (1, "V"), "time": None, "period": 60.0},
+            ],
+        )
+
     def test_read_strings_repeat(self):
         experiment = pybamm.Experiment(
             ["Discharge at 10 mA for 0.5 hours"]
@@ -120,11 +173,11 @@ def test_read_strings_repeat(self):
         self.assertEqual(
             experiment.operating_conditions,
             [
-                (0.01, "A", 1800.0, 60),
-                (-0.5, "C", 2700.0, 60),
-                (1, "V", 20.0, 60),
-                (-0.5, "C", 2700.0, 60),
-                (1, "V", 20.0, 60),
+                {"electric": (0.01, "A"), "time": 1800.0, "period": 60},
+                {"electric": (-0.5, "C"), "time": 2700.0, "period": 60},
+                {"electric": (1, "V"), "time": 20.0, "period": 60},
+                {"electric": (-0.5, "C"), "time": 2700.0, "period": 60},
+                {"electric": (1, "V"), "time": 20.0, "period": 60},
             ],
         )
         self.assertEqual(experiment.period, 60)
@@ -140,10 +193,10 @@ def test_cycle_unpacking(self):
         self.assertEqual(
             experiment.operating_conditions,
             [
-                (0.05, "C", 1800.0, 60.0),
-                (-0.2, "C", 2700.0, 60.0),
-                (0.05, "C", 1800.0, 60.0),
-                (-0.2, "C", 2700.0, 60.0),
+                {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0},
+                {"electric": (-0.2, "C"), "time": 2700.0, "period": 60.0},
+                {"electric": (0.05, "C"), "time": 1800.0, "period": 60.0},
+                {"electric": (-0.2, "C"), "time": 2700.0, "period": 60.0},
             ],
         )
         self.assertEqual(experiment.cycle_lengths, [2, 1, 1])
@@ -159,6 +212,11 @@ def test_str_repr(self):
         )
 
     def test_bad_strings(self):
+        with self.assertRaisesRegex(ValueError, "cccv_handling"):
+            pybamm.Experiment(
+                ["Discharge at 1 C for 20 seconds", "Charge at 0.5 W for 10 minutes"],
+                cccv_handling="bad",
+            )
         with self.assertRaisesRegex(
             TypeError, "Operating conditions should be strings or tuples of strings"
         ):
diff --git a/tests/unit/test_experiments/test_simulation_with_experiment.py b/tests/unit/test_experiments/test_simulation_with_experiment.py
index 32e932c071..2f5d11d859 100644
--- a/tests/unit/test_experiments/test_simulation_with_experiment.py
+++ b/tests/unit/test_experiments/test_simulation_with_experiment.py
@@ -117,6 +117,47 @@ def test_run_experiment(self):
         self.assertEqual(len(sol3.cycles), 2)
         os.remove("test_experiment.sav")
 
+    def test_run_experiment_cccv_ode(self):
+        experiment_2step = pybamm.Experiment(
+            [
+                (
+                    "Discharge at C/20 for 1 hour",
+                    "Charge at 1 A until 4.1 V",
+                    "Hold at 4.1 V until C/2",
+                    "Discharge at 2 W for 1 hour",
+                ),
+            ],
+        )
+        experiment_ode = pybamm.Experiment(
+            [
+                (
+                    "Discharge at C/20 for 1 hour",
+                    "Charge at 1 A until 4.1 V",
+                    "Hold at 4.1 V until C/2",
+                    "Discharge at 2 W for 1 hour",
+                ),
+            ],
+            cccv_handling="ode",
+        )
+        solutions = []
+        for experiment in [experiment_2step, experiment_ode]:
+            model = pybamm.lithium_ion.SPM()
+            sim = pybamm.Simulation(model, experiment=experiment)
+            solution = sim.solve(solver=pybamm.CasadiSolver("fast with events"))
+            solutions.append(solution)
+
+        np.testing.assert_array_almost_equal(
+            solutions[0]["Terminal voltage [V]"].data,
+            solutions[1]["Terminal voltage [V]"].data,
+            decimal=2,
+        )
+        np.testing.assert_array_almost_equal(
+            solutions[0]["Current [A]"].data,
+            solutions[1]["Current [A]"].data,
+            decimal=0,
+        )
+        self.assertEqual(solutions[1].termination, "final time")
+
     def test_run_experiment_old_setup_type(self):
         experiment = pybamm.Experiment(
             [
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index a00c5a91ba..c77082de92 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -25,7 +25,7 @@
 'lithium plating': 'none' (possible: ['none', 'reversible', 'irreversible'])
 'lithium plating porosity change': 'false' (possible: ['true', 'false'])
 'loss of active material': 'stress-driven' (possible: ['none', 'stress-driven', 'reaction-driven'])
-'operating mode': 'current' (possible: ['current', 'voltage', 'power'])
+'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'CCCV'])
 'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'fast diffusion', 'uniform profile', 'quadratic profile', 'quartic profile'])
 'particle mechanics': 'swelling only' (possible: ['none', 'swelling only', 'swelling and cracking'])
 'particle shape': 'spherical' (possible: ['spherical', 'user', 'no particles'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
index f2dd194eb8..0e5960f8a1 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_spm.py
@@ -191,6 +191,11 @@ def test_well_posed_power(self):
         model = pybamm.lithium_ion.SPM(options)
         model.check_well_posedness()
 
+    def test_well_posed_cccv(self):
+        options = {"operating mode": "CCCV"}
+        model = pybamm.lithium_ion.SPM(options)
+        model.check_well_posedness()
+
     def test_well_posed_function(self):
         def external_circuit_function(variables):
             I = variables["Current [A]"]
diff --git a/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py b/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py
index 9cfd51cf99..9aa743030e 100644
--- a/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py
+++ b/tests/unit/test_models/test_submodels/test_external_circuit/test_function_control.py
@@ -31,6 +31,13 @@ def test_public_functions(self):
         std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
+    def test_cccv_control(self):
+        param = pybamm.LithiumIonParameters()
+        submodel = pybamm.external_circuit.CCCVFunctionControl(param)
+        variables = {"Terminal voltage [V]": pybamm.Scalar(0)}
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")