diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 75d7f1427..993d4d37e 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -393,7 +393,10 @@ class Simulation:
             # to correspond to a single discharge
             elif t_eval is None:
                 C_rate = self._parameter_values["C-rate"]
-                t_end = 3600 / C_rate
+                try:
+                    t_end = 3600 / C_rate
+                except TypeError:
+                    t_end = 3600
                 t_eval = np.linspace(0, t_end, 100)
 
             self.t_eval = t_eval
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 50099b7e4..8dd2ce273 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -219,21 +219,31 @@ class BaseSolver(object):
 
         # Check for heaviside functions in rhs and algebraic and add discontinuity
         # events if these exist.
-        # Note: only checks for the case of t < X, t <= X, X < t, or X <= t
-        for symbol in itertools.chain(
-            model.concatenated_rhs.pre_order(), model.concatenated_algebraic.pre_order()
-        ):
-            if isinstance(symbol, pybamm.Heaviside):
-                if symbol.right.id == pybamm.t.id:
-                    expr = symbol.left
-                elif symbol.left.id == pybamm.t.id:
-                    expr = symbol.right
-
-                model.events.append(
-                    pybamm.Event(
-                        str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY
+        # Note: only checks for the case of t < X, t <= X, X < t, or X <= t, but also
+        # accounts for the fact that t might be dimensional
+        # Only do this for DAE models as ODE models can deal with discontinuities fine
+        if len(model.algebraic) > 0:
+            for symbol in itertools.chain(
+                model.concatenated_rhs.pre_order(),
+                model.concatenated_algebraic.pre_order(),
+            ):
+                if isinstance(symbol, pybamm.Heaviside):
+                    # Dimensionless
+                    if symbol.right.id == pybamm.t.id:
+                        expr = symbol.left
+                    elif symbol.left.id == pybamm.t.id:
+                        expr = symbol.right
+                    # Dimensional
+                    elif symbol.right.id == (pybamm.t * model.timescale).id:
+                        expr = symbol.left.new_copy() / symbol.right.right.new_copy()
+                    elif symbol.left.id == (pybamm.t * model.timescale).id:
+                        expr = symbol.right.new_copy() / symbol.left.right.new_copy()
+
+                    model.events.append(
+                        pybamm.Event(
+                            str(symbol), expr.new_copy(), pybamm.EventType.DISCONTINUITY
+                        )
                     )
-                )
 
         # Process rhs, algebraic and event expressions
         rhs, rhs_eval, jac_rhs = process(model.concatenated_rhs, "RHS")
@@ -328,7 +338,7 @@ class BaseSolver(object):
             Initial conditions that are consistent with the algebraic equations (roots
             of the algebraic equations)
         """
-        pybamm.logger.info("Start calculating consistent initial conditions")
+        pybamm.logger.info("Start calculating consistent states")
         rhs = model.rhs_eval
         algebraic = model.algebraic_eval
         jac = model.jac_algebraic_eval
@@ -336,7 +346,7 @@ class BaseSolver(object):
             y0_guess = model.concatenated_initial_conditions.flatten()
 
         # Split y0_guess into differential and algebraic
-        len_rhs = rhs(0, y0_guess).shape[0]
+        len_rhs = rhs(time, y0_guess).shape[0]
         y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs])
 
         def root_fun(y0_alg):
@@ -344,8 +354,8 @@ class BaseSolver(object):
             y0 = np.concatenate([y0_diff, y0_alg])
             out = algebraic(time, y0)
             pybamm.logger.debug(
-                "Evaluating algebraic equations at t=0, L2-norm is {}".format(
-                    np.linalg.norm(out)
+                "Evaluating algebraic equations at t={}, L2-norm is {}".format(
+                    time * model.timescale, np.linalg.norm(out)
                 )
             )
             return out
@@ -472,43 +482,46 @@ class BaseSolver(object):
 
         # make sure they are increasing in time
         discontinuities = sorted(discontinuities)
-        if len(discontinuities) > 0:
-            pybamm.logger.info(
-                "Discontinuity events found at t = {}".format(discontinuities)
-            )
-        else:
-            pybamm.logger.info("No discontinuity events found")
 
         # remove any identical discontinuities
         discontinuities = [
             v
             for i, v in enumerate(discontinuities)
-            if i == len(discontinuities) - 1
-            or discontinuities[i] < discontinuities[i + 1]
+            if (
+                i == len(discontinuities) - 1
+                or discontinuities[i] < discontinuities[i + 1]
+            )
+            and v > 0
         ]
 
+        if len(discontinuities) > 0:
+            pybamm.logger.info(
+                "Discontinuity events found at t = {}".format(discontinuities)
+            )
+        else:
+            pybamm.logger.info("No discontinuity events found")
+
         # insert time points around discontinuities in t_eval
         # keep track of sub sections to integrate by storing start and end indices
         start_indices = [0]
         end_indices = []
+        eps = sys.float_info.epsilon
         for dtime in discontinuities:
             dindex = np.searchsorted(t_eval_dimensionless, dtime, side="left")
             end_indices.append(dindex + 1)
             start_indices.append(dindex + 1)
-            if t_eval_dimensionless[dindex] == dtime:
-                t_eval_dimensionless[dindex] += sys.float_info.epsilon
+            if dtime - eps < t_eval_dimensionless[dindex] < dtime + eps:
+                t_eval_dimensionless[dindex] += eps
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless, dindex, dtime - sys.float_info.epsilon
+                    t_eval_dimensionless, dindex, dtime - eps
                 )
             else:
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless,
-                    dindex,
-                    [dtime - sys.float_info.epsilon, dtime + sys.float_info.epsilon],
+                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps],
                 )
         end_indices.append(len(t_eval_dimensionless))
 
-        # integrate separatly over each time segment and accumulate into the solution
+        # integrate separately over each time segment and accumulate into the solution
         # object, restarting the solver at each discontinuity (and recalculating a
         # consistent state afterwards if a dae)
         old_y0 = model.y0
@@ -521,16 +534,13 @@ class BaseSolver(object):
                 )
             )
             timer.reset()
+            new_solution = self._integrate(
+                model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
+            )
+            new_solution.solve_time = timer.time()
             if solution is None:
-                solution = self._integrate(
-                    model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
-                )
-                solution.solve_time = timer.time()
+                solution = new_solution
             else:
-                new_solution = self._integrate(
-                    model, t_eval_dimensionless[start_index:end_index], ext_and_inputs
-                )
-                new_solution.solve_time = timer.time()
                 solution.append(new_solution, start_index=0)
 
             if solution.termination != "final time":
@@ -538,14 +548,6 @@ class BaseSolver(object):
 
             if end_index != len(t_eval_dimensionless):
                 # setup for next integration subsection
-                y0_guess = solution.y[:, -1]
-                if model.algebraic:
-                    model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], y0_guess
-                    )
-                else:
-                    model.y0 = y0_guess
-
                 last_state = solution.y[:, -1]
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 141364ac2..ea809fda7 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -100,13 +100,13 @@ class CasadiSolver(pybamm.BaseSolver):
             return solution
         elif self.mode == "safe":
             # Step-and-check
+            t = t_eval[0]
             init_event_signs = np.sign(
                 np.concatenate(
-                    [event(0, model.y0) for event in model.terminate_events_eval]
+                    [event(t, model.y0) for event in model.terminate_events_eval]
                 )
             )
             pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
-            t = t_eval[0]
             y0 = model.y0
             # Initialize solution
             solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis])
@@ -144,7 +144,7 @@ class CasadiSolver(pybamm.BaseSolver):
                 new_event_signs = np.sign(
                     np.concatenate(
                         [
-                            event(0, current_step_sol.y[:, -1])
+                            event(t, current_step_sol.y[:, -1])
                             for event in model.terminate_events_eval
                         ]
                     )
@@ -186,15 +186,15 @@ class CasadiSolver(pybamm.BaseSolver):
             # set up and solve
             t = casadi.MX.sym("t")
             u = casadi.MX.sym("u", u_stacked.shape[0])
-            y_diff = casadi.MX.sym("y_diff", rhs(0, y0, u).shape[0])
+            y_diff = casadi.MX.sym("y_diff", rhs(t_eval[0], y0, u).shape[0])
             problem = {"t": t, "x": y_diff, "p": u}
-            if algebraic(0, y0, u).is_empty():
+            if algebraic(t_eval[0], y0, u).is_empty():
                 method = "cvodes"
                 problem.update({"ode": rhs(t, y_diff, u)})
             else:
                 options["calc_ic"] = True
                 method = "idas"
-                y_alg = casadi.MX.sym("y_alg", algebraic(0, y0, u).shape[0])
+                y_alg = casadi.MX.sym("y_alg", algebraic(t_eval[0], y0, u).shape[0])
                 y_full = casadi.vertcat(y_diff, y_alg)
                 problem.update(
                     {
@@ -215,7 +215,7 @@ class CasadiSolver(pybamm.BaseSolver):
         )
 
     def _run_integrator(self, integrator, model, y0, inputs, t_eval):
-        rhs_size = model.rhs_eval(0, y0).shape[0]
+        rhs_size = model.rhs_eval(t_eval[0], y0).shape[0]
         y0_diff, y0_alg = np.split(y0, [rhs_size])
         try:
             # Try solving

Problem with simulating step changes in current
Car current simulation doesn't work. Seems to not pick up the right discontinuities because time is dimensional
diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py
index af9c6a179..0d87c632e 100644
--- a/tests/unit/test_simulation.py
+++ b/tests/unit/test_simulation.py
@@ -385,6 +385,29 @@ class TestSimulation(unittest.TestCase):
         with self.assertWarns(pybamm.SolverWarning):
             sim.solve(t_eval=np.linspace(0, time_data[-1], 800))
 
+    def test_discontinuous_current(self):
+        def car_current(t):
+            current = (
+                1 * (t >= 0) * (t <= 1000)
+                - 0.5 * (1000 < t) * (t <= 2000)
+                + 0.5 * (2000 < t)
+            )
+            return current
+
+        model = pybamm.lithium_ion.DFN()
+        param = model.default_parameter_values
+        param["Current function [A]"] = car_current
+
+        sim = pybamm.Simulation(
+            model, parameter_values=param, solver=pybamm.CasadiSolver(mode="fast")
+        )
+        sim.solve()
+        current = sim.solution["Current [A]"]
+        tau = model.timescale.evaluate()
+        self.assertEqual(current(0), 1)
+        self.assertEqual(current(1500 / tau), -0.5)
+        self.assertEqual(current(3000 / tau), 0.5)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 3921f2031..62120bcde 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -50,6 +50,7 @@ class TestBaseSolver(unittest.TestCase):
         class ScalarModel:
             concatenated_initial_conditions = np.array([[2]])
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -68,6 +69,7 @@ class TestBaseSolver(unittest.TestCase):
         class VectorModel:
             concatenated_initial_conditions = np.zeros_like(vec)
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return y[0:1]
@@ -101,6 +103,7 @@ class TestBaseSolver(unittest.TestCase):
         class Model:
             concatenated_initial_conditions = np.array([2])
             jac_algebraic_eval = None
+            timescale = 1
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -123,6 +126,16 @@ class TestBaseSolver(unittest.TestCase):
         ):
             solver.calculate_consistent_state(Model())
 
+    def test_time_too_short(self):
+        solver = pybamm.BaseSolver()
+        model = pybamm.BaseModel()
+        v = pybamm.StateVector(slice(0, 1))
+        model.rhs = {v: v}
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "It looks like t_eval might be dimensionless"
+        ):
+            solver.solve(model, np.linspace(0, 0.1))
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 7e6974c57..da2e6a3a6 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -259,6 +259,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         def nonsmooth_mult(t):
             return int(t < discontinuity) + 1.0
+
         rate = pybamm.Function(nonsmooth_rate, pybamm.t)
         mult = pybamm.Function(nonsmooth_mult, pybamm.t)
         model.rhs = {var1: rate * var1}
@@ -267,14 +268,16 @@ class TestScikitsSolvers(unittest.TestCase):
         model.events = [
             pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
             pybamm.Event("var2 = 2.5", pybamm.min(var2 - 2.5)),
-            pybamm.Event("nonsmooth rate",
-                         pybamm.Scalar(discontinuity),
-                         pybamm.EventType.DISCONTINUITY
-                         ),
-            pybamm.Event("nonsmooth mult",
-                         pybamm.Scalar(discontinuity),
-                         pybamm.EventType.DISCONTINUITY
-                         )
+            pybamm.Event(
+                "nonsmooth rate",
+                pybamm.Scalar(discontinuity),
+                pybamm.EventType.DISCONTINUITY,
+            ),
+            pybamm.Event(
+                "nonsmooth mult",
+                pybamm.Scalar(discontinuity),
+                pybamm.EventType.DISCONTINUITY,
+            ),
         ]
         disc = get_discretisation_for_testing()
         disc.process_model(model)
@@ -285,9 +288,9 @@ class TestScikitsSolvers(unittest.TestCase):
         # create two time series, one without a time point on the discontinuity,
         # and one with
         t_eval1 = np.linspace(0, 5, 10)
-        t_eval2 = np.insert(t_eval1,
-                            np.searchsorted(t_eval1, discontinuity),
-                            discontinuity)
+        t_eval2 = np.insert(
+            t_eval1, np.searchsorted(t_eval1, discontinuity), discontinuity
+        )
         solution1 = solver.solve(model, t_eval1)
         solution2 = solver.solve(model, t_eval2)
 
@@ -312,13 +315,13 @@ class TestScikitsSolvers(unittest.TestCase):
             np.testing.assert_array_less(solution.y[-1], 2.5)
             var1_soln = np.exp(0.2 * solution.t)
             y0 = np.exp(0.2 * discontinuity)
-            var1_soln[solution.t > discontinuity] = \
-                y0 * np.exp(
+            var1_soln[solution.t > discontinuity] = y0 * np.exp(
                 0.1 * (solution.t[solution.t > discontinuity] - discontinuity)
             )
             var2_soln = 2 * var1_soln
-            var2_soln[solution.t > discontinuity] = \
-                var1_soln[solution.t > discontinuity]
+            var2_soln[solution.t > discontinuity] = var1_soln[
+                solution.t > discontinuity
+            ]
             np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
             np.testing.assert_allclose(solution.y[-1], var2_soln, rtol=1e-06)
 
@@ -592,6 +595,95 @@ class TestScikitsSolvers(unittest.TestCase):
             step_solution.y[-1], 2 * np.exp(0.1 * step_solution.t), decimal=5
         )
 
+    def test_model_solver_dae_nonsmooth(self):
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2")
+        discontinuity = 0.6
+
+        # Create three different models with the same solution, each expressing the
+        # discontinuity in a different way
+
+        # first model explicitly adds a discontinuity event
+        def nonsmooth_rate(t):
+            return 0.1 * (t < discontinuity) + 0.1
+
+        rate = pybamm.Function(nonsmooth_rate, pybamm.t)
+        model1 = pybamm.BaseModel()
+        model1.rhs = {var1: rate * var1}
+        model1.algebraic = {var2: var2}
+        model1.initial_conditions = {var1: 1, var2: 0}
+        model1.events = [
+            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
+            pybamm.Event(
+                "nonsmooth rate",
+                pybamm.Scalar(discontinuity),
+                pybamm.EventType.DISCONTINUITY,
+            ),
+        ]
+
+        # second model implicitly adds a discontinuity event via a heaviside function
+        model2 = pybamm.BaseModel()
+        model2.rhs = {var1: (0.1 * (pybamm.t < discontinuity) + 0.1) * var1}
+        model2.algebraic = {var2: var2}
+        model2.initial_conditions = {var1: 1, var2: 0}
+        model2.events = [
+            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
+        ]
+
+        # third model implicitly adds a discontinuity event via another heaviside
+        # function
+        model3 = pybamm.BaseModel()
+        model3.rhs = {var1: (-0.1 * (discontinuity < pybamm.t) + 0.2) * var1}
+        model3.algebraic = {var2: var2}
+        model3.initial_conditions = {var1: 1, var2: 0}
+        model3.events = [
+            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
+        ]
+
+        for model in [model1, model2, model3]:
+
+            disc = get_discretisation_for_testing()
+            disc.process_model(model)
+
+            # Solve
+            solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+
+            # create two time series, one without a time point on the discontinuity,
+            # and one with
+            t_eval1 = np.linspace(0, 5, 10)
+            t_eval2 = np.insert(
+                t_eval1, np.searchsorted(t_eval1, discontinuity), discontinuity
+            )
+            solution1 = solver.solve(model, t_eval1)
+            solution2 = solver.solve(model, t_eval2)
+
+            # check time vectors
+            for solution in [solution1, solution2]:
+                # time vectors are ordered
+                self.assertTrue(np.all(solution.t[:-1] <= solution.t[1:]))
+
+                # time value before and after discontinuity is an epsilon away
+                dindex = np.searchsorted(solution.t, discontinuity)
+                value_before = solution.t[dindex - 1]
+                value_after = solution.t[dindex]
+                self.assertEqual(value_before + sys.float_info.epsilon, discontinuity)
+                self.assertEqual(value_after - sys.float_info.epsilon, discontinuity)
+
+            # both solution time vectors should have same number of points
+            self.assertEqual(len(solution1.t), len(solution2.t))
+
+            # check solution
+            for solution in [solution1, solution2]:
+                np.testing.assert_array_less(solution.y[0], 1.5)
+                np.testing.assert_array_less(solution.y[-1], 2.5)
+                var1_soln = np.exp(0.2 * solution.t)
+                y0 = np.exp(0.2 * discontinuity)
+                var1_soln[solution.t > discontinuity] = y0 * np.exp(
+                    0.1 * (solution.t[solution.t > discontinuity] - discontinuity)
+                )
+                np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
+
     def test_ode_solver_fail_with_dae(self):
         model = pybamm.BaseModel()
         a = pybamm.Scalar(1)
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index dbd2aabc0..fe064b64c 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -6,7 +6,6 @@ import unittest
 import numpy as np
 from tests import get_mesh_for_testing
 import warnings
-from tests import get_discretisation_for_testing
 import sys
 
 
@@ -138,91 +137,6 @@ class TestScipySolver(unittest.TestCase):
             np.ones((N, T.size)) * (T[np.newaxis, :] - np.exp(T[np.newaxis, :])),
         )
 
-    def test_model_solver_ode_nonsmooth(self):
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        var1 = pybamm.Variable("var1", domain=whole_cell)
-        discontinuity = 0.6
-
-        # Create three different models with the same solution, each expressing the
-        # discontinuity in a different way
-
-        # first model explicitly adds a discontinuity event
-        def nonsmooth_rate(t):
-            return 0.1 * (t < discontinuity) + 0.1
-
-        rate = pybamm.Function(nonsmooth_rate, pybamm.t)
-        model1 = pybamm.BaseModel()
-        model1.rhs = {var1: rate * var1}
-        model1.initial_conditions = {var1: 1}
-        model1.events = [
-            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
-            pybamm.Event(
-                "nonsmooth rate",
-                pybamm.Scalar(discontinuity),
-                pybamm.EventType.DISCONTINUITY,
-            ),
-        ]
-
-        # second model implicitly adds a discontinuity event via a heaviside function
-        model2 = pybamm.BaseModel()
-        model2.rhs = {var1: (0.1 * (pybamm.t < discontinuity) + 0.1) * var1}
-        model2.initial_conditions = {var1: 1}
-        model2.events = [
-            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
-        ]
-
-        # third model implicitly adds a discontinuity event via another heaviside
-        # function
-        model3 = pybamm.BaseModel()
-        model3.rhs = {var1: (-0.1 * (discontinuity < pybamm.t) + 0.2) * var1}
-        model3.initial_conditions = {var1: 1}
-        model3.events = [
-            pybamm.Event("var1 = 1.5", pybamm.min(var1 - 1.5)),
-        ]
-
-        for model in [model1, model2, model3]:
-
-            disc = get_discretisation_for_testing()
-            disc.process_model(model)
-
-            # Solve
-            solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8)
-
-            # create two time series, one without a time point on the discontinuity,
-            # and one with
-            t_eval1 = np.linspace(0, 5, 10)
-            t_eval2 = np.insert(
-                t_eval1, np.searchsorted(t_eval1, discontinuity), discontinuity
-            )
-            solution1 = solver.solve(model, t_eval1)
-            solution2 = solver.solve(model, t_eval2)
-
-            # check time vectors
-            for solution in [solution1, solution2]:
-                # time vectors are ordered
-                self.assertTrue(np.all(solution.t[:-1] <= solution.t[1:]))
-
-                # time value before and after discontinuity is an epsilon away
-                dindex = np.searchsorted(solution.t, discontinuity)
-                value_before = solution.t[dindex - 1]
-                value_after = solution.t[dindex]
-                self.assertEqual(value_before + sys.float_info.epsilon, discontinuity)
-                self.assertEqual(value_after - sys.float_info.epsilon, discontinuity)
-
-            # both solution time vectors should have same number of points
-            self.assertEqual(len(solution1.t), len(solution2.t))
-
-            # check solution
-            for solution in [solution1, solution2]:
-                np.testing.assert_array_less(solution.y[0], 1.5)
-                np.testing.assert_array_less(solution.y[-1], 2.5)
-                var1_soln = np.exp(0.2 * solution.t)
-                y0 = np.exp(0.2 * discontinuity)
-                var1_soln[solution.t > discontinuity] = y0 * np.exp(
-                    0.1 * (solution.t[solution.t > discontinuity] - discontinuity)
-                )
-                np.testing.assert_allclose(solution.y[0], var1_soln, rtol=1e-06)
-
     def test_model_step_python(self):
         # Create model
         model = pybamm.BaseModel()

pip install -e .[dev]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gcc gfortran
3.9
requirements-docs.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
anyio==4.9.0
anytree==2.12.1
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
asttokens==3.0.0
async-lru==2.0.5
attrs==25.3.0
autograd==1.7.0
babel==2.17.0
beautifulsoup4==4.13.3
black==25.1.0
bleach==6.2.0
casadi==3.7.0
certifi==2025.1.31
cffi==1.17.1
charset-normalizer==3.4.1
click==8.1.8
comm==0.2.2
contourpy==1.3.0
cycler==0.12.1
debugpy==1.8.13
decorator==5.2.1
defusedxml==0.7.1
exceptiongroup==1.2.2
executing==2.2.0
fastjsonschema==2.21.1
flake8==7.2.0
fonttools==4.56.0
fqdn==1.5.1
h11==0.14.0
httpcore==1.0.7
httpx==0.28.1
idna==3.10
importlib_metadata==8.6.1
importlib_resources==6.5.2
iniconfig==2.1.0
ipykernel==6.29.5
ipython==8.18.1
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.10.0
jsonpointer==3.0.0
jsonschema==4.23.0
jsonschema-specifications==2024.10.1
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.12.0
jupyter-lsp==2.2.5
jupyter_client==8.6.3
jupyter_core==5.7.2
jupyter_server==2.15.0
jupyter_server_terminals==0.5.3
jupyterlab==4.3.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.27.3
jupyterlab_widgets==3.0.13
kiwisolver==1.4.7
Levenshtein==0.27.1
MarkupSafe==3.0.2
matplotlib==3.9.4
matplotlib-inline==0.1.7
mccabe==0.7.0
mistune==3.1.3
mypy-extensions==1.0.0
nbclient==0.10.2
nbconvert==7.16.6
nbformat==5.10.4
nest-asyncio==1.6.0
notebook==7.3.3
notebook_shim==0.2.4
numpy==2.0.2
overrides==7.7.0
packaging==24.2
pandas==2.2.3
pandocfilters==1.5.1
parso==0.8.4
pathspec==0.12.1
pexpect==4.9.0
pillow==11.1.0
platformdirs==4.3.7
pluggy==1.5.0
prometheus_client==0.21.1
prompt_toolkit==3.0.50
psutil==7.0.0
ptyprocess==0.7.0
pure_eval==0.2.3
-e git+https://github.com/pybamm-team/PyBaMM.git@084d8e4ca3bef12aa11e5d7546de6817f9bd3a9c#egg=pybamm
pycodestyle==2.13.0
pycparser==2.22
pyflakes==3.3.1
Pygments==2.19.1
pyparsing==3.2.3
pytest==8.3.5
python-dateutil==2.9.0.post0
python-json-logger==3.3.0
python-Levenshtein==0.27.1
pytz==2025.2
PyYAML==6.0.2
pyzmq==26.3.0
RapidFuzz==3.12.2
referencing==0.36.2
requests==2.32.3
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.24.0
scikit-fem==10.0.2
scipy==1.13.1
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
soupsieve==2.6
stack-data==0.6.3
terminado==0.18.1
tinycss2==1.4.0
tomli==2.2.1
tornado==6.4.2
traitlets==5.14.3
types-python-dateutil==2.9.0.20241206
typing_extensions==4.13.0
tzdata==2025.2
uri-template==1.3.0
urllib3==2.3.0
wcwidth==0.2.13
webcolors==24.11.1
webencodings==0.5.1
websocket-client==1.8.0
widgetsnbextension==4.0.13
zipp==3.21.0

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

tests/unit/test_simulation.py::TestSimulation::test_discontinuous_current
tests/unit/test_simulation.py::TestSimulation::test_basic_ops tests/unit/test_simulation.py::TestSimulation::test_drive_cycle_data tests/unit/test_simulation.py::TestSimulation::test_get_variable_array tests/unit/test_simulation.py::TestSimulation::test_plot tests/unit/test_simulation.py::TestSimulation::test_reuse_commands tests/unit/test_simulation.py::TestSimulation::test_save_load tests/unit/test_simulation.py::TestSimulation::test_save_load_dae tests/unit/test_simulation.py::TestSimulation::test_set_crate tests/unit/test_simulation.py::TestSimulation::test_set_defaults tests/unit/test_simulation.py::TestSimulation::test_set_defaults2 tests/unit/test_simulation.py::TestSimulation::test_set_external_variable tests/unit/test_simulation.py::TestSimulation::test_solve tests/unit/test_simulation.py::TestSimulation::test_specs tests/unit/test_simulation.py::TestSimulation::test_step tests/unit/test_simulation.py::TestSimulation::test_step_with_inputs tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_base_solver_init tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_fail_consistent_initial_conditions tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_find_consistent_initial_conditions tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_nonmonotonic_teval tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_ode_solver_fail_with_dae tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_step_or_solve_empty_model tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_time_too_short tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_failure tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_ode_with_jacobian_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_event_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_event_with_casadi tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_external tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_inputs tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_inputs_with_casadi tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_step_python
pybamm/solvers/casadi_solver.py pybamm/simulation.py pybamm/solvers/base_solver.py
pybamm/solvers/casadi_solver.py pybamm/simulation.py pybamm/solvers/base_solver.py