diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3e6cf962d..6f2e22249 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844))
 -   Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824))
 -   Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Added options to export a solution to matlab or csv ([#811](https://github.com/pybamm-team/PyBaMM/pull/811))
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
index fd01abfa2..16d17eef2 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
@@ -177,17 +177,11 @@ class BasicDFN(BaseModel):
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T),
-                "Neumann",
-            ),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T),
-                "Neumann",
-            ),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
         }
         # c_n_init and c_p_init can in general be functions of x
         # Note the broadcasting, for domains
@@ -201,22 +195,14 @@ class BasicDFN(BaseModel):
         self.initial_conditions[c_s_p] = param.c_p_init(x_p)
         # Events specify points at which a solution should terminate
         self.events += [
-            pybamm.Event(
-                "Minimum negative particle surface concentration",
-                pybamm.min(c_s_surf_n) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum negative particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_n),
-            ),
-            pybamm.Event(
-                "Minimum positive particle surface concentration",
-                pybamm.min(c_s_surf_p) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum positive particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_p),
-            ),
+            pybamm.Event("Minimum negative particle surface concentration",
+                         pybamm.min(c_s_surf_n) - 0.01),
+            pybamm.Event("Maximum negative particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
+            pybamm.Event("Minimum positive particle surface concentration",
+                         pybamm.min(c_s_surf_p) - 0.01),
+            pybamm.Event("Maximum positive particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
         ]
         ######################
         # Current in the solid
@@ -271,11 +257,8 @@ class BasicDFN(BaseModel):
             "right": (pybamm.Scalar(0), "Neumann"),
         }
         self.initial_conditions[c_e] = param.c_e_init
-        self.events.append(
-            pybamm.Event(
-                "Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002
-            )
-        )
+        self.events.append(pybamm.Event("Zero electrolyte concentration cut-off",
+                                        pybamm.min(c_e) - 0.002))
 
         ######################
         # (Some) variables
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
index 60956d210..15d8e2ab9 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
@@ -80,48 +80,34 @@ class BasicSPM(BaseModel):
         N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p)
         self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n)
         self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p)
-        # Surf takes the surface value of a variable, i.e. its boundary value on the
-        # right side. This is also accessible via `boundary_value(x, "right")`, with
-        # "left" providing the boundary value of the left side
-        c_s_surf_n = pybamm.surf(c_s_n)
-        c_s_surf_p = pybamm.surf(c_s_p)
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T),
-                "Neumann",
-            ),
+            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (
-                -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T),
-                "Neumann",
-            ),
+            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
         }
         # c_n_init and c_p_init are functions, but for the SPM we evaluate them at x=0
         # and x=1 since there is no x-dependence in the particles
         self.initial_conditions[c_s_n] = param.c_n_init(0)
         self.initial_conditions[c_s_p] = param.c_p_init(1)
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        c_s_surf_p = pybamm.surf(c_s_p)
         # Events specify points at which a solution should terminate
         self.events += [
-            pybamm.Event(
-                "Minimum negative particle surface concentration",
-                pybamm.min(c_s_surf_n) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum negative particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_n),
-            ),
-            pybamm.Event(
-                "Minimum positive particle surface concentration",
-                pybamm.min(c_s_surf_p) - 0.01,
-            ),
-            pybamm.Event(
-                "Maximum positive particle surface concentration",
-                (1 - 0.01) - pybamm.max(c_s_surf_p),
-            ),
+            pybamm.Event("Minimum negative particle surface concentration",
+                         pybamm.min(c_s_surf_n) - 0.01),
+            pybamm.Event("Maximum negative particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
+            pybamm.Event("Minimum positive particle surface concentration",
+                         pybamm.min(c_s_surf_p) - 0.01),
+            pybamm.Event("Maximum positive particle surface concentration",
+                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
         ]
 
         # Note that the SPM does not have any algebraic equations, so the `algebraic`
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 1c031b8b8..806015e90 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -23,7 +23,10 @@ class BaseSolver(object):
     atol : float, optional
         The absolute tolerance for the solver (default is 1e-6).
     root_method : str, optional
-        The method to use to find initial conditions (default is "lm")
+        The method to use to find initial conditions (default is "casadi"). If "casadi",
+        the solver uses casadi's Newton rootfinding algorithm to find initial
+        conditions. Otherwise, the solver uses 'scipy.optimize.root' with method
+        specified by 'root_method' (e.g. "lm", "hybr", ...)
     root_tol : float, optional
         The tolerance for the initial-condition solver (default is 1e-6).
     max_steps: int, optional
@@ -36,7 +39,7 @@ class BaseSolver(object):
         method=None,
         rtol=1e-6,
         atol=1e-6,
-        root_method="lm",
+        root_method="casadi",
         root_tol=1e-6,
         max_steps=1000,
     ):
@@ -125,10 +128,11 @@ class BaseSolver(object):
                 "Cannot use ODE solver '{}' to solve DAE model".format(self.name)
             )
 
+        if self.ode_solver is True:
+            self.root_method = None
         if (
-            isinstance(self, pybamm.CasadiSolver)
-            and model.convert_to_format != "casadi"
-        ):
+            isinstance(self, pybamm.CasadiSolver) or self.root_method == "casadi"
+        ) and model.convert_to_format != "casadi":
             pybamm.logger.warning(
                 f"Converting {model.name} to CasADi for solving with CasADi solver"
             )
@@ -268,6 +272,16 @@ class BaseSolver(object):
         model.terminate_events_eval = terminate_events_eval
         model.discontinuity_events_eval = discontinuity_events_eval
 
+        # Save CasADi functions for the CasADi solver
+        # Note: when we pass to casadi the ode part of the problem must be in explicit
+        # form so we pre-multiply by the inverse of the mass matrix
+        if self.root_method == "casadi" or isinstance(self, pybamm.CasadiSolver):
+            mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries)
+            explicit_rhs = mass_matrix_inv @ rhs(t_casadi, y_casadi, u_casadi_stacked)
+            model.casadi_rhs = casadi.Function(
+                "rhs", [t_casadi, y_casadi, u_casadi_stacked], [explicit_rhs]
+            )
+            model.casadi_algebraic = algebraic
         # Calculate consistent initial conditions for the algebraic equations
         if len(model.algebraic) > 0:
             all_states = pybamm.NumpyConcatenation(
@@ -278,26 +292,13 @@ class BaseSolver(object):
             model.residuals_eval = residuals_eval
             model.jacobian_eval = jacobian_eval
             y0_guess = y0.flatten()
-            model.y0 = self.calculate_consistent_state(model, 0, y0_guess)
+            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
         else:
             # can use DAE solver to solve ODE model
             model.residuals_eval = Residuals(rhs, "residuals", model)
             model.jacobian_eval = jac_rhs
             model.y0 = y0.flatten()
 
-        # Save CasADi functions for the CasADi solver
-        # Note: when we pass to casadi the ode part of the problem must be in explicit
-        # form so we pre-multiply by the inverse of the mass matrix
-        if model.convert_to_format == "casadi" and isinstance(
-            self, pybamm.CasadiSolver
-        ):
-            mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries)
-            explicit_rhs = mass_matrix_inv @ rhs(t_casadi, y_casadi, u_casadi_stacked)
-            model.casadi_rhs = casadi.Function(
-                "rhs", [t_casadi, y_casadi, u_casadi_stacked], [explicit_rhs]
-            )
-            model.casadi_algebraic = algebraic
-
         pybamm.logger.info("Finish solver set-up")
 
     def set_inputs(self, model, ext_and_inputs):
@@ -319,7 +320,7 @@ class BaseSolver(object):
         if model.jacobian_eval:
             model.jacobian_eval.set_inputs(ext_and_inputs)
 
-    def calculate_consistent_state(self, model, time=0, y0_guess=None):
+    def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None):
         """
         Calculate consistent state for the algebraic equations through
         root-finding
@@ -328,6 +329,12 @@ class BaseSolver(object):
         ----------
         model : :class:`pybamm.BaseModel`
             The model for which to calculate initial conditions.
+        time : float
+            The time at which to calculate the states
+        y0_guess : :class:`np.array`
+            Guess for the rootfinding
+        inputs : dict, optional
+            Any input parameters to pass to the model when solving
 
         Returns
         -------
@@ -336,65 +343,101 @@ class BaseSolver(object):
             of the algebraic equations)
         """
         pybamm.logger.info("Start calculating consistent states")
-        rhs = model.rhs_eval
-        algebraic = model.algebraic_eval
-        jac = model.jac_algebraic_eval
         if y0_guess is None:
             y0_guess = model.concatenated_initial_conditions.flatten()
 
         # Split y0_guess into differential and algebraic
-        len_rhs = rhs(time, y0_guess).shape[0]
+        len_rhs = model.rhs_eval(time, y0_guess).shape[0]
         y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs])
+        inputs = inputs or {}
 
-        def root_fun(y0_alg):
-            "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)"
-            y0 = np.concatenate([y0_diff, y0_alg])
-            out = algebraic(time, y0)
-            pybamm.logger.debug(
-                "Evaluating algebraic equations at t={}, L2-norm is {}".format(
-                    time * model.timescale, np.linalg.norm(out)
-                )
+        # Solve using casadi or scipy
+        if self.root_method == "casadi":
+            # Set up
+            u_stacked = casadi.vertcat(*[x for x in inputs.values()])
+            u = casadi.MX.sym("u", u_stacked.shape[0])
+            y_alg = casadi.MX.sym("y_alg", y0_alg_guess.shape[0])
+            y = casadi.vertcat(y0_diff, y_alg)
+            alg_root = model.casadi_algebraic(time, y, u)
+            # Solve
+            # set error_on_fail to False and just check the final output is small
+            # enough
+            roots = casadi.rootfinder(
+                "roots",
+                "newton",
+                dict(x=y_alg, p=u, g=alg_root),
+                {"abstol": self.root_tol},
             )
-            return out
-
-        if jac:
-            if issparse(jac(0, y0_guess)):
+            try:
+                y0_alg = roots(y0_alg_guess, u_stacked).full().flatten()
+                success = True
+                message = None
+                # Check final output
+                fun = model.casadi_algebraic(
+                    time, casadi.vertcat(y0_diff, y0_alg), u_stacked
+                )
+            except RuntimeError as err:
+                success = False
+                message = err.args[0]
+                fun = None
+        else:
+            algebraic = model.algebraic_eval
+            jac = model.jac_algebraic_eval
+
+            def root_fun(y0_alg):
+                "Evaluates algebraic using y0_diff (fixed) and y0_alg (changed by algo)"
+                y0 = np.concatenate([y0_diff, y0_alg])
+                out = algebraic(time, y0)
+                pybamm.logger.debug(
+                    "Evaluating algebraic equations at t={}, L2-norm is {}".format(
+                        time * model.timescale, np.linalg.norm(out)
+                    )
+                )
+                return out
 
-                def jac_fn(y0_alg):
-                    """
-                    Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
-                    """
-                    y0 = np.concatenate([y0_diff, y0_alg])
-                    return jac(0, y0)[:, len_rhs:].toarray()
+            if jac:
+                if issparse(jac(0, y0_guess)):
 
-            else:
+                    def jac_fn(y0_alg):
+                        """
+                        Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
+                        """
+                        y0 = np.concatenate([y0_diff, y0_alg])
+                        return jac(0, y0)[:, len_rhs:].toarray()
 
-                def jac_fn(y0_alg):
-                    """
-                    Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
-                    """
-                    y0 = np.concatenate([y0_diff, y0_alg])
-                    return jac(0, y0)[:, len_rhs:]
+                else:
 
-        else:
-            jac_fn = None
-        # Find the values of y0_alg that are roots of the algebraic equations
-        sol = optimize.root(
-            root_fun,
-            y0_alg_guess,
-            jac=jac_fn,
-            method=self.root_method,
-            tol=self.root_tol,
-        )
-        # Return full set of consistent initial conditions (y0_diff unchanged)
-        y0_consistent = np.concatenate([y0_diff, sol.x])
+                    def jac_fn(y0_alg):
+                        """
+                        Evaluates jacobian using y0_diff (fixed) and y0_alg (varying)
+                        """
+                        y0 = np.concatenate([y0_diff, y0_alg])
+                        return jac(0, y0)[:, len_rhs:]
 
-        if sol.success and np.all(sol.fun < self.root_tol * len(sol.x)):
+            else:
+                jac_fn = None
+            # Find the values of y0_alg that are roots of the algebraic equations
+            sol = optimize.root(
+                root_fun,
+                y0_alg_guess,
+                jac=jac_fn,
+                method=self.root_method,
+                tol=self.root_tol,
+            )
+            # Set outputs
+            y0_alg = sol.x
+            success = sol.success
+            fun = sol.fun
+            message = sol.message
+
+        if success and np.all(fun < self.root_tol * len(y0_alg)):
+            # Return full set of consistent initial conditions (y0_diff unchanged)
+            y0_consistent = np.concatenate([y0_diff, y0_alg])
             pybamm.logger.info("Finish calculating consistent initial conditions")
             return y0_consistent
-        elif not sol.success:
+        elif not success:
             raise pybamm.SolverError(
-                "Could not find consistent initial conditions: {}".format(sol.message)
+                "Could not find consistent initial conditions: {}".format(message)
             )
         else:
             raise pybamm.SolverError(
@@ -402,7 +445,7 @@ class BaseSolver(object):
                 Could not find consistent initial conditions: solver terminated
                 successfully, but maximum solution error ({}) above tolerance ({})
                 """.format(
-                    np.max(sol.fun), self.root_tol * len(sol.x)
+                    np.max(fun), self.root_tol * len(y0_alg)
                 )
             )
 
@@ -514,7 +557,7 @@ class BaseSolver(object):
                 )
             else:
                 t_eval_dimensionless = np.insert(
-                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps],
+                    t_eval_dimensionless, dindex, [dtime - eps, dtime + eps]
                 )
         end_indices.append(len(t_eval_dimensionless))
 
@@ -548,7 +591,10 @@ class BaseSolver(object):
                 last_state = solution.y[:, -1]
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
-                        model, t_eval_dimensionless[end_index], last_state
+                        model,
+                        t_eval_dimensionless[end_index],
+                        last_state,
+                        ext_and_inputs,
                     )
                 else:
                     model.y0 = last_state
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index ea809fda7..34fd74043 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -47,7 +47,7 @@ class CasadiSolver(pybamm.BaseSolver):
         mode="safe",
         rtol=1e-6,
         atol=1e-6,
-        root_method="lm",
+        root_method="casadi",
         root_tol=1e-6,
         max_step_decrease_count=5,
         **extra_options,
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 795f3361b..be4a1932c 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -36,7 +36,7 @@ class IDAKLUSolver(pybamm.BaseSolver):
     """
 
     def __init__(
-        self, rtol=1e-6, atol=1e-6, root_method="lm", root_tol=1e-6, max_steps=1000
+        self, rtol=1e-6, atol=1e-6, root_method="casadi", root_tol=1e-6, max_steps=1000
     ):
 
         if idaklu_spec is None:
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index 5f9999cbb..74620cafe 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -40,7 +40,7 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
         method="ida",
         rtol=1e-6,
         atol=1e-6,
-        root_method="lm",
+        root_method="casadi",
         root_tol=1e-6,
         max_steps=1000,
     ):
diff --git a/setup.py b/setup.py
index db141d4f4..fa13bd203 100644
--- a/setup.py
+++ b/setup.py
@@ -492,7 +492,7 @@ setup(
     # List of dependencies
     install_requires=[
         "numpy>=1.16",
-        "scipy>=1.0",
+        "scipy>=1.3",
         "pandas>=0.24",
         "anytree>=2.4.3",
         "autograd>=1.2",

Initialize with casadi
**Summary**
Use casadi root-finder to initialize DAEs

 **Motivation**
 Should be faster



diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index ce44197a3..7ab105982 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -1,6 +1,7 @@
 #
 # Tests for the Base Solver class
 #
+import casadi
 import pybamm
 import numpy as np
 from scipy.sparse import csr_matrix
@@ -49,9 +50,16 @@ class TestBaseSolver(unittest.TestCase):
     def test_find_consistent_initial_conditions(self):
         # Simple system: a single algebraic equation
         class ScalarModel:
-            concatenated_initial_conditions = np.array([[2]])
-            jac_algebraic_eval = None
-            timescale = 1
+            def __init__(self):
+                self.concatenated_initial_conditions = np.array([[2]])
+                self.jac_algebraic_eval = None
+                self.timescale = 1
+                t = casadi.MX.sym("t")
+                y = casadi.MX.sym("y")
+                u = casadi.MX.sym("u")
+                self.casadi_algebraic = casadi.Function(
+                    "alg", [t, y, u], [self.algebraic_eval(t, y)]
+                )
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -59,18 +67,30 @@ class TestBaseSolver(unittest.TestCase):
             def algebraic_eval(self, t, y):
                 return y + 2
 
-        solver = pybamm.BaseSolver()
+        solver = pybamm.BaseSolver(root_method="lm")
         model = ScalarModel()
         init_cond = solver.calculate_consistent_state(model)
         np.testing.assert_array_equal(init_cond, -2)
+        # with casadi
+        solver_with_casadi = pybamm.BaseSolver(root_method="casadi", root_tol=1e-12)
+        model = ScalarModel()
+        init_cond = solver_with_casadi.calculate_consistent_state(model)
+        np.testing.assert_array_equal(init_cond, -2)
 
         # More complicated system
         vec = np.array([0.0, 1.0, 1.5, 2.0])
 
         class VectorModel:
-            concatenated_initial_conditions = np.zeros_like(vec)
-            jac_algebraic_eval = None
-            timescale = 1
+            def __init__(self):
+                self.concatenated_initial_conditions = np.zeros_like(vec)
+                self.jac_algebraic_eval = None
+                self.timescale = 1
+                t = casadi.MX.sym("t")
+                y = casadi.MX.sym("y", vec.size)
+                u = casadi.MX.sym("u")
+                self.casadi_algebraic = casadi.Function(
+                    "alg", [t, y, u], [self.algebraic_eval(t, y)]
+                )
 
             def rhs_eval(self, t, y):
                 return y[0:1]
@@ -81,6 +101,9 @@ class TestBaseSolver(unittest.TestCase):
         model = VectorModel()
         init_cond = solver.calculate_consistent_state(model)
         np.testing.assert_array_almost_equal(init_cond, vec)
+        # with casadi
+        init_cond = solver_with_casadi.calculate_consistent_state(model)
+        np.testing.assert_array_almost_equal(init_cond, vec)
 
         # With jacobian
         def jac_dense(t, y):
@@ -102,9 +125,16 @@ class TestBaseSolver(unittest.TestCase):
 
     def test_fail_consistent_initial_conditions(self):
         class Model:
-            concatenated_initial_conditions = np.array([2])
-            jac_algebraic_eval = None
-            timescale = 1
+            def __init__(self):
+                self.concatenated_initial_conditions = np.array([2])
+                self.jac_algebraic_eval = None
+                self.timescale = 1
+                t = casadi.MX.sym("t")
+                y = casadi.MX.sym("y")
+                u = casadi.MX.sym("u")
+                self.casadi_algebraic = casadi.Function(
+                    "alg", [t, y, u], [self.algebraic_eval(t, y)]
+                )
 
             def rhs_eval(self, t, y):
                 return np.array([])
@@ -120,12 +150,19 @@ class TestBaseSolver(unittest.TestCase):
             "Could not find consistent initial conditions: The iteration is not making",
         ):
             solver.calculate_consistent_state(Model())
-        solver = pybamm.BaseSolver()
+        solver = pybamm.BaseSolver(root_method="lm")
         with self.assertRaisesRegex(
             pybamm.SolverError,
             "Could not find consistent initial conditions: solver terminated",
         ):
             solver.calculate_consistent_state(Model())
+        # with casadi
+        solver = pybamm.BaseSolver(root_method="casadi")
+        with self.assertRaisesRegex(
+            pybamm.SolverError,
+            "Could not find consistent initial conditions: .../casadi",
+        ):
+            solver.calculate_consistent_state(Model())
 
     def test_time_too_short(self):
         solver = pybamm.BaseSolver()
diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py
index 9683b5cbc..703673ae8 100644
--- a/tests/unit/test_solvers/test_idaklu_solver.py
+++ b/tests/unit/test_solvers/test_idaklu_solver.py
@@ -28,7 +28,7 @@ class TestIDAKLUSolver(unittest.TestCase):
             disc = pybamm.Discretisation()
             disc.process_model(model)
 
-            solver = pybamm.IDAKLUSolver()
+            solver = pybamm.IDAKLUSolver(root_method="lm")
 
             t_eval = np.linspace(0, 3, 100)
             solution = solver.solve(model, t_eval)
@@ -58,7 +58,7 @@ class TestIDAKLUSolver(unittest.TestCase):
         mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
         disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
         disc.process_model(model)
-        solver = pybamm.IDAKLUSolver()
+        solver = pybamm.IDAKLUSolver(root_method="lm")
 
         variable_tols = {"Electrolyte concentration": 1e-3}
         solver.set_atol_by_variable(variable_tols, model)
@@ -76,7 +76,7 @@ class TestIDAKLUSolver(unittest.TestCase):
         disc = pybamm.Discretisation()
         disc.process_model(model)
 
-        solver = pybamm.IDAKLUSolver()
+        solver = pybamm.IDAKLUSolver(root_method="lm")
 
         t_eval = np.linspace(0, 3, 100)
         with self.assertRaisesRegex(pybamm.SolverError, "KLU requires the Jacobian"):
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index e3b39791a..905507c69 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -194,7 +194,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -214,7 +214,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -238,7 +238,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 5, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_less(solution.y[0], 1.5)
@@ -283,7 +283,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
 
         # create two time series, one without a time point on the discontinuity,
         # and one with
@@ -355,7 +355,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         model.jacobian = jacobian
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -372,7 +372,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc.process_model(model)
 
         # Solve
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -421,7 +421,7 @@ class TestScikitsSolvers(unittest.TestCase):
         disc = get_discretisation_for_testing()
         disc.process_model(model)
 
-        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
 
         # Step once
         dt = 1
@@ -514,7 +514,10 @@ class TestScikitsSolvers(unittest.TestCase):
             disc.process_model(model)
 
             # Solve
-            solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+            if form == "python":
+                solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8, root_method="lm")
+            else:
+                solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
             t_eval = np.linspace(0, 5, 100)
             solution = solver.solve(model, t_eval, inputs={"rate 1": 0.1, "rate 2": 2})
             np.testing.assert_array_less(solution.y[0], 1.5)

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y wget gcc gfortran libopenblas-dev
3.7
requirements/base.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.7.1
anytree==2.12.1
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
attrs==24.2.0
autograd==1.6.2
Babel==2.14.0
backcall==0.2.0
beautifulsoup4==4.13.3
black==23.3.0
bleach==6.0.0
casadi==3.7.0
certifi @ file:///croot/certifi_1671487769961/work/certifi
cffi==1.15.1
charset-normalizer==3.4.1
click==8.1.8
comm==0.1.4
cycler==0.11.0
debugpy==1.7.0
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
exceptiongroup==1.2.2
fastjsonschema==2.21.1
flake8==5.0.4
fonttools==4.38.0
future==1.0.0
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
importlib-metadata==4.2.0
importlib-resources==5.12.0
iniconfig==2.0.0
ipykernel==6.16.2
ipython==7.34.0
ipython-genutils==0.2.0
ipywidgets==8.1.5
jedi==0.19.2
Jinja2==3.1.6
json5==0.9.16
jsonschema==4.17.3
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-server==1.24.0
jupyter_client==7.4.9
jupyter_core==4.12.0
jupyterlab==3.5.3
jupyterlab-pygments==0.2.2
jupyterlab_server==2.15.2
jupyterlab_widgets==3.0.13
kiwisolver==1.4.5
Levenshtein==0.23.0
MarkupSafe==2.1.5
matplotlib==3.5.3
matplotlib-inline==0.1.6
mccabe==0.7.0
mistune==3.0.2
mypy-extensions==1.0.0
nbclassic==1.2.0
nbclient==0.7.4
nbconvert==7.6.0
nbformat==5.8.0
nest-asyncio==1.6.0
notebook==6.5.7
notebook_shim==0.2.4
numpy==1.21.6
packaging==24.0
pandas==1.3.5
pandocfilters==1.5.1
parso==0.8.4
pathspec==0.11.2
pexpect==4.9.0
pickleshare==0.7.5
Pillow==9.5.0
pkgutil_resolve_name==1.3.10
platformdirs==4.0.0
pluggy==1.2.0
prometheus-client==0.17.1
prompt_toolkit==3.0.48
psutil==7.0.0
ptyprocess==0.7.0
-e git+https://github.com/pybamm-team/PyBaMM.git@11080ac2373ffb1b378ac65355eb300607334ff7#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.17.2
pyparsing==3.1.4
pyrsistent==0.19.3
pytest==7.4.4
python-dateutil==2.9.0.post0
python-Levenshtein==0.23.0
pytz==2025.2
pyzmq==26.2.1
rapidfuzz==3.4.0
requests==2.31.0
scikit-fem==8.1.0
scipy==1.7.3
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.4.1
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.17.1
tinycss2==1.2.1
tomli==2.0.1
tornado==6.2
traitlets==5.9.0
typed-ast==1.5.5
typing_extensions==4.7.1
urllib3==2.0.7
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.6.1
widgetsnbextension==4.0.13
zipp==3.15.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2022.12.7=py37h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=22.3.1=py37h06a4308_0
  - python=3.7.16=h7a1cb2a_0
  - readline=8.2=h5eee18b_0
  - setuptools=65.6.3=py37h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.38.4=py37h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.7.1
      - anytree==2.12.1
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - attrs==24.2.0
      - autograd==1.6.2
      - babel==2.14.0
      - backcall==0.2.0
      - beautifulsoup4==4.13.3
      - black==23.3.0
      - bleach==6.0.0
      - casadi==3.7.0
      - cffi==1.15.1
      - charset-normalizer==3.4.1
      - click==8.1.8
      - comm==0.1.4
      - cycler==0.11.0
      - debugpy==1.7.0
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - exceptiongroup==1.2.2
      - fastjsonschema==2.21.1
      - flake8==5.0.4
      - fonttools==4.38.0
      - future==1.0.0
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==4.2.0
      - importlib-resources==5.12.0
      - iniconfig==2.0.0
      - ipykernel==6.16.2
      - ipython==7.34.0
      - ipython-genutils==0.2.0
      - ipywidgets==8.1.5
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.9.16
      - jsonschema==4.17.3
      - jupyter==1.1.1
      - jupyter-client==7.4.9
      - jupyter-console==6.6.3
      - jupyter-core==4.12.0
      - jupyter-server==1.24.0
      - jupyterlab==3.5.3
      - jupyterlab-pygments==0.2.2
      - jupyterlab-server==2.15.2
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.5
      - levenshtein==0.23.0
      - markupsafe==2.1.5
      - matplotlib==3.5.3
      - matplotlib-inline==0.1.6
      - mccabe==0.7.0
      - mistune==3.0.2
      - mypy-extensions==1.0.0
      - nbclassic==1.2.0
      - nbclient==0.7.4
      - nbconvert==7.6.0
      - nbformat==5.8.0
      - nest-asyncio==1.6.0
      - notebook==6.5.7
      - notebook-shim==0.2.4
      - numpy==1.21.6
      - packaging==24.0
      - pandas==1.3.5
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pathspec==0.11.2
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==9.5.0
      - pkgutil-resolve-name==1.3.10
      - platformdirs==4.0.0
      - pluggy==1.2.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.48
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.17.2
      - pyparsing==3.1.4
      - pyrsistent==0.19.3
      - pytest==7.4.4
      - python-dateutil==2.9.0.post0
      - python-levenshtein==0.23.0
      - pytz==2025.2
      - pyzmq==26.2.1
      - rapidfuzz==3.4.0
      - requests==2.31.0
      - scikit-fem==8.1.0
      - scipy==1.7.3
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - snowballstemmer==2.2.0
      - soupsieve==2.4.1
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.17.1
      - tinycss2==1.2.1
      - tomli==2.0.1
      - tornado==6.2
      - traitlets==5.9.0
      - typed-ast==1.5.5
      - typing-extensions==4.7.1
      - urllib3==2.0.7
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.6.1
      - widgetsnbextension==4.0.13
      - zipp==3.15.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_find_consistent_initial_conditions
tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_fail_consistent_initial_conditions
tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_base_solver_init 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
pybamm/solvers/casadi_solver.py pybamm/models/full_battery_models/lithium_ion/basic_spm.py setup.py CHANGELOG.md pybamm/solvers/idaklu_solver.py pybamm/solvers/base_solver.py pybamm/solvers/scikits_dae_solver.py pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
pybamm/solvers/casadi_solver.py pybamm/models/full_battery_models/lithium_ion/basic_spm.py setup.py CHANGELOG.md pybamm/solvers/idaklu_solver.py pybamm/solvers/base_solver.py pybamm/solvers/scikits_dae_solver.py pybamm/models/full_battery_models/lithium_ion/basic_dfn.py