diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py
index 21ae70e8c..da93e50d9 100644
--- a/pybamm/solvers/algebraic_solver.py
+++ b/pybamm/solvers/algebraic_solver.py
@@ -5,11 +5,13 @@ import casadi
 import pybamm
 import numpy as np
 from scipy import optimize
+from scipy.sparse import issparse
 
 
 class AlgebraicSolver(pybamm.BaseSolver):
     """Solve a discretised model which contains only (time independent) algebraic
     equations using a root finding algorithm.
+    Uses scipy.optimize.root.
     Note: this solver could be extended for quasi-static models, or models in
     which the time derivative is manually discretised and results in a (possibly
     nonlinear) algebaric system at each time level.
@@ -20,11 +22,16 @@ class AlgebraicSolver(pybamm.BaseSolver):
         The method to use to solve the system (default is "lm")
     tol : float, optional
         The tolerance for the solver (default is 1e-6).
+    extra_options : dict, optional
+        Any options to pass to the rootfinder. Vary depending on which method is chosen.
+        Please consult `SciPy documentation <https://tinyurl.com/ybr6cfqs>`_ for
+        details.
     """
 
-    def __init__(self, method="lm", tol=1e-6):
+    def __init__(self, method="lm", tol=1e-6, extra_options=None):
         super().__init__(method=method)
         self.tol = tol
+        self.extra_options = extra_options or {}
         self.name = "Algebraic solver ({})".format(method)
         self.algebraic_solver = True
         pybamm.citations.register("virtanen2020scipy")
@@ -54,47 +61,82 @@ class AlgebraicSolver(pybamm.BaseSolver):
         if model.convert_to_format == "casadi":
             inputs = casadi.vertcat(*[x for x in inputs.values()])
 
-        algebraic = model.algebraic_eval
         y0 = model.y0
+        if isinstance(y0, casadi.DM):
+            y0 = y0.full().flatten()
+
+        # The casadi algebraic solver can read rhs equations, but leaves them unchanged
+        # i.e. the part of the solution vector that corresponds to the differential
+        # equations will be equal to the initial condition provided. This allows this
+        # solver to be used for initialising the DAE solvers
+        # Split y0 into differential and algebraic
+        if model.rhs == {}:
+            len_rhs = 0
+        else:
+            len_rhs = model.rhs_eval(t_eval[0], y0, inputs).shape[0]
+        y0_diff, y0_alg = np.split(y0, [len_rhs])
+
+        algebraic = model.algebraic_eval
 
-        y = np.empty((len(y0), len(t_eval)))
+        y_alg = np.empty((len(y0_alg), len(t_eval)))
 
         for idx, t in enumerate(t_eval):
 
-            def root_fun(y):
+            def root_fun(y_alg):
                 "Evaluates algebraic using y"
+                y = np.concatenate([y0_diff, y_alg])
                 out = algebraic(t, y, inputs)
                 pybamm.logger.debug(
                     "Evaluating algebraic equations at t={}, L2-norm is {}".format(
-                        t, np.linalg.norm(out)
+                        t * model.timescale_eval, np.linalg.norm(out)
                     )
                 )
                 return out
 
-            if model.jacobian_eval is not None:
+            jac = model.jac_algebraic_eval
+            if jac:
+                if issparse(jac(t_eval[0], y0, inputs)):
+
+                    def jac_fn(y_alg):
+                        """
+                        Evaluates jacobian using y0_diff (fixed) and y_alg (varying)
+                        """
+                        y = np.concatenate([y0_diff, y_alg])
+                        return jac(0, y, inputs)[:, len_rhs:].toarray()
+
+                else:
 
-                def jac(y):
-                    return model.jacobian_eval(t, y, inputs)
+                    def jac_fn(y_alg):
+                        """
+                        Evaluates jacobian using y0_diff (fixed) and y_alg (varying)
+                        """
+                        y = np.concatenate([y0_diff, y_alg])
+                        return jac(0, y, inputs)[:, len_rhs:]
 
             else:
-                jac = None
+                jac_fn = None
 
             # Evaluate algebraic with new t and previous y0, if it's already close
             # enough then keep it
             if np.all(abs(algebraic(t, y0, inputs)) < self.tol):
                 pybamm.logger.debug("Keeping same solution at t={}".format(t))
-                y[:, idx] = y0
+                y_alg[:, idx] = y0_alg
             # Otherwise calculate new y0
             else:
                 sol = optimize.root(
-                    root_fun, y0, method=self.method, tol=self.tol, jac=jac
+                    root_fun,
+                    y0_alg,
+                    method=self.method,
+                    tol=self.tol,
+                    jac=jac_fn,
+                    options=self.extra_options,
                 )
 
                 if sol.success and np.all(abs(sol.fun) < self.tol):
                     # update initial guess for the next iteration
-                    y0 = sol.x
+                    y0_alg = sol.x
                     # update solution array
-                    y[:, idx] = y0
+                    y_alg[:, idx] = y0_alg
                 elif not sol.success:
                     raise pybamm.SolverError(
                         "Could not find acceptable solution: {}".format(sol.message)
@@ -110,5 +152,8 @@ class AlgebraicSolver(pybamm.BaseSolver):
                         )
                     )
 
+        # Concatenate differential part
+        y_diff = np.r_[[y0_diff] * len(t_eval)].T
+        y_sol = np.r_[y_diff, y_alg]
         # Return solution object (no events, so pass None to t_event, y_event)
-        return pybamm.Solution(t_eval, y, termination="success")
+        return pybamm.Solution(t_eval, y_sol, termination="success")
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index d76e66dca..6f4db433d 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -6,8 +6,6 @@ import copy
 import pybamm
 import numbers
 import numpy as np
-from scipy import optimize
-from scipy.sparse import issparse
 import sys
 import itertools
 
@@ -23,16 +21,15 @@ class BaseSolver(object):
         The relative tolerance for the solver (default is 1e-6).
     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 "casadi"). If "casadi",
+    root_method : str or pybamm algebraic solver class, optional
+        The method to use to find initial conditions (for DAE solvers).
+        If a solver class, must be an algebraic solver class.
+        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
-        The maximum number of steps the solver will take before terminating
-        (default is 1000).
     """
 
     def __init__(
@@ -40,17 +37,20 @@ class BaseSolver(object):
         method=None,
         rtol=1e-6,
         atol=1e-6,
-        root_method="casadi",
+        root_method=None,
         root_tol=1e-6,
-        max_steps=1000,
+        max_steps="deprecated",
     ):
         self._method = method
         self._rtol = rtol
         self._atol = atol
-        self.root_method = root_method
         self.root_tol = root_tol
-        self.max_steps = max_steps
-
+        self.root_method = root_method
+        if max_steps != "deprecated":
+            raise ValueError(
+                "max_steps has been deprecated, and should be set using the "
+                "solver-specific extra-options dictionaries instead"
+            )
         self.models_set_up = set()
 
         # Defaults, can be overwritten by specific solver
@@ -88,6 +88,18 @@ class BaseSolver(object):
 
     @root_method.setter
     def root_method(self, method):
+        if method == "casadi":
+            method = pybamm.CasadiAlgebraicSolver(self.root_tol)
+        elif isinstance(method, str):
+            method = pybamm.AlgebraicSolver(method, self.root_tol)
+        elif not (
+            method is None
+            or (
+                isinstance(method, pybamm.BaseSolver)
+                and method.algebraic_solver is True
+            )
+        ):
+            raise pybamm.SolverError("Root method must be an algebraic solver")
         self._root_method = method
 
     @property
@@ -98,14 +110,6 @@ class BaseSolver(object):
     def root_tol(self, tol):
         self._root_tol = tol
 
-    @property
-    def max_steps(self):
-        return self._max_steps
-
-    @max_steps.setter
-    def max_steps(self, max_steps):
-        self._max_steps = max_steps
-
     def copy(self):
         "Returns a copy of the solver"
         new_solver = copy.copy(self)
@@ -136,6 +140,13 @@ class BaseSolver(object):
             raise pybamm.SolverError(
                 """Cannot use algebraic solver to solve model with time derivatives"""
             )
+        # casadi solver won't allow solving algebraic model so we have to raise an
+        # error here
+        if isinstance(self, pybamm.CasadiSolver) and len(model.rhs) == 0:
+            raise pybamm.SolverError(
+                "Cannot use CasadiSolver to solve algebraic model, "
+                "use CasadiAlgebraicSolver instead"
+            )
         # Discretise model if it isn't already discretised
         # This only works with purely 0D models, as otherwise the mesh and spatial
         # method should be specified by the user
@@ -150,15 +161,10 @@ class BaseSolver(object):
                 )
 
         inputs = inputs or {}
-        model.y0 = model.concatenated_initial_conditions.evaluate(
-            0, None, inputs=inputs
-        ).flatten()
 
         # Set model timescale
         model.timescale_eval = model.timescale.evaluate(inputs=inputs)
 
-        if self.ode_solver is True:
-            self.root_method = None
         if (
             isinstance(self, (pybamm.CasadiSolver, pybamm.CasadiAlgebraicSolver))
         ) and model.convert_to_format != "casadi":
@@ -168,7 +174,10 @@ class BaseSolver(object):
                 )
             )
             model.convert_to_format = "casadi"
-        if self.root_method == "casadi" and model.convert_to_format != "casadi":
+        if (
+            isinstance(self.root_method, pybamm.CasadiAlgebraicSolver)
+            and model.convert_to_format != "casadi"
+        ):
             pybamm.logger.warning(
                 "Converting {} to CasADi for calculating ICs with CasADi".format(
                     model.name
@@ -179,35 +188,14 @@ class BaseSolver(object):
         if model.convert_to_format != "casadi":
             simp = pybamm.Simplification()
             # Create Jacobian from concatenated rhs and algebraic
-            y = pybamm.StateVector(slice(0, np.size(model.y0)))
+            y = pybamm.StateVector(slice(0, model.concatenated_initial_conditions.size))
             # set up Jacobian object, for re-use of dict
             jacobian = pybamm.Jacobian()
         else:
-            # Create placeholder inputs for evaluating rhs and algebraic sizes
-            placeholder_inputs = {}
-            for k, v in inputs.items():
-                if isinstance(v, casadi.MX):
-                    placeholder_inputs[k] = np.zeros(v.shape[0])
-                else:
-                    placeholder_inputs[k] = v
             # Convert model attributes to casadi
             t_casadi = casadi.MX.sym("t")
-            y_diff = casadi.MX.sym(
-                "y_diff",
-                len(
-                    model.concatenated_rhs.evaluate(
-                        0, model.y0, inputs=placeholder_inputs
-                    )
-                ),
-            )
-            y_alg = casadi.MX.sym(
-                "y_alg",
-                len(
-                    model.concatenated_algebraic.evaluate(
-                        0, model.y0, inputs=placeholder_inputs
-                    )
-                ),
-            )
+            y_diff = casadi.MX.sym("y_diff", model.concatenated_rhs.size)
+            y_alg = casadi.MX.sym("y_alg", model.concatenated_algebraic.size)
             y_casadi = casadi.vertcat(y_diff, y_alg)
             p_casadi = {}
             for name, value in inputs.items():
@@ -308,6 +296,14 @@ class BaseSolver(object):
                             )
                         )
 
+        # Process initial conditions
+        initial_conditions = process(
+            model.concatenated_initial_conditions,
+            "initial_conditions",
+            use_jacobian=False,
+        )[0]
+        init_eval = InitialConditions(initial_conditions, model)
+
         # Process rhs, algebraic and event expressions
         rhs, rhs_eval, jac_rhs = process(model.concatenated_rhs, "RHS")
         algebraic, algebraic_eval, jac_algebraic = process(
@@ -328,16 +324,22 @@ class BaseSolver(object):
         ]
 
         # Add the solver attributes
+        model.init_eval = init_eval
         model.rhs_eval = rhs_eval
         model.algebraic_eval = algebraic_eval
         model.jac_algebraic_eval = jac_algebraic
         model.terminate_events_eval = terminate_events_eval
         model.discontinuity_events_eval = discontinuity_events_eval
 
+        # Calculate initial conditions
+        model.y0 = init_eval(inputs)
+
         # 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):
+        if isinstance(self.root_method, pybamm.CasadiAlgebraicSolver) or isinstance(
+            self, (pybamm.CasadiSolver, pybamm.CasadiAlgebraicSolver)
+        ):
             # can use DAE solver to solve model with algebraic equations only
             if len(model.rhs) > 0:
                 mass_matrix_inv = casadi.MX(model.mass_matrix_inv.entries)
@@ -386,35 +388,35 @@ class BaseSolver(object):
 
         """
         if self.algebraic_solver is True:
+            # Don't update model.y0
             return None
         elif len(model.algebraic) == 0:
             if update_rhs is True:
                 # Recalculate initial conditions for the rhs equations
-                model.y0 = model.concatenated_initial_conditions.evaluate(
-                    0, None, inputs=inputs
-                ).flatten()
+                model.y0 = model.init_eval(inputs)
             else:
+                # Don't update model.y0
                 return None
         else:
             if update_rhs is True:
                 # Recalculate initial conditions for the rhs equations
-                y0_from_inputs = model.concatenated_initial_conditions.evaluate(
-                    0, None, inputs=inputs
-                ).flatten()
+                y0_from_inputs = model.init_eval(inputs)
                 # Reuse old solution for algebraic equations
                 y0_from_model = model.y0
-                len_rhs = len(
-                    model.concatenated_rhs.evaluate(0, model.y0, inputs=inputs)
-                )
-                y0_guess = np.r_[y0_from_inputs[:len_rhs], y0_from_model[len_rhs:]]
-            else:
-                y0_guess = model.y0
-            model.y0 = self.calculate_consistent_state(model, 0, y0_guess, inputs)
+                len_rhs = model.concatenated_rhs.size
+                # update model.y0, which is used for initialising the algebraic solver
+                if len_rhs == 0:
+                    model.y0 = y0_from_model
+                else:
+                    model.y0 = casadi.vertcat(
+                        y0_from_inputs[:len_rhs], y0_from_model[len_rhs:]
+                    )
+            model.y0 = self.calculate_consistent_state(model, 0, inputs)
 
-    def calculate_consistent_state(self, model, time=0, y0_guess=None, inputs=None):
+    def calculate_consistent_state(self, model, time=0, inputs=None):
         """
         Calculate consistent state for the algebraic equations through
-        root-finding
+        root-finding. model.y0 is used as the initial guess for rootfinding
 
         Parameters
         ----------
@@ -422,8 +424,6 @@ class BaseSolver(object):
             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
 
@@ -434,115 +434,13 @@ class BaseSolver(object):
             of the algebraic equations)
         """
         pybamm.logger.info("Start calculating consistent states")
-        if y0_guess is None:
-            y0_guess = model.concatenated_initial_conditions.flatten()
-
-        inputs = inputs or {}
-        if model.convert_to_format == "casadi":
-            inputs = casadi.vertcat(*[x for x in inputs.values()])
-
-        # Split y0_guess into differential and algebraic
-        len_rhs = model.rhs_eval(time, y0_guess, inputs).shape[0]
-        y0_diff, y0_alg_guess = np.split(y0_guess, [len_rhs])
-
-        # Solve using casadi or scipy
-        if self.root_method == "casadi":
-            # Set up
-            p = casadi.MX.sym("p", inputs.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, p)
-            # Solve
-            roots = casadi.rootfinder(
-                "roots",
-                "newton",
-                dict(x=y_alg, p=p, g=alg_root),
-                {"abstol": self.root_tol},
-            )
-            try:
-                y0_alg = roots(y0_alg_guess, inputs).full().flatten()
-                success = True
-                message = None
-                # Check final output
-                fun = model.casadi_algebraic(
-                    time, casadi.vertcat(y0_diff, y0_alg), inputs
-                )
-                abs_fun = casadi.fabs(fun)
-                max_fun = casadi.mmax(fun)
-            except RuntimeError as err:
-                success = False
-                message = err.args[0]
-                abs_fun = None
-                max_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, inputs)
-                pybamm.logger.debug(
-                    "Evaluating algebraic equations at t={}, L2-norm is {}".format(
-                        time * model.timescale, np.linalg.norm(out)
-                    )
-                )
-                return out
-
-            if jac:
-                if issparse(jac(0, y0_guess, inputs)):
-
-                    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, inputs)[:, len_rhs:].toarray()
-
-                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, inputs)[:, len_rhs:]
-
-            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,
-            )
-            pybamm.citations.register("virtanen2020scipy")
-
-            # Set outputs
-            y0_alg = sol.x
-            success = sol.success
-            fun = sol.fun
-            abs_fun = np.abs(fun)
-            max_fun = np.max(fun)
-            message = sol.message
-
-        if success and np.all(abs_fun < self.root_tol):
-            # 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 success:
-            raise pybamm.SolverError(
-                "Could not find consistent initial conditions: {}".format(message)
-            )
-        else:
+        try:
+            root_sol = self.root_method._integrate(model, [time], inputs)
+        except pybamm.SolverError as e:
             raise pybamm.SolverError(
-                "Could not find consistent initial conditions: solver terminated "
-                "successfully, but maximum solution error "
-                "({}) above tolerance ({})".format(max_fun, self.root_tol)
+                "Could not find consistent initial conditions: {}".format(e.args[0])
             )
+        return root_sol.y.flatten()
 
     def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         """
@@ -693,15 +591,13 @@ class BaseSolver(object):
             if end_index != len(t_eval_dimensionless):
                 # setup for next integration subsection
                 last_state = solution.y[:, -1]
+                # update y0 (for DAE solvers, this updates the initial guess for the
+                # rootfinder)
+                model.y0 = last_state
                 if len(model.algebraic) > 0:
                     model.y0 = self.calculate_consistent_state(
-                        model,
-                        t_eval_dimensionless[end_index],
-                        last_state,
-                        ext_and_inputs,
+                        model, t_eval_dimensionless[end_index], ext_and_inputs,
                     )
-                else:
-                    model.y0 = last_state
 
         # restore old y0
         model.y0 = old_y0
@@ -954,3 +850,19 @@ class Residuals(SolverCallable):
     def __call__(self, t, y, ydot, inputs):
         states_eval = super().__call__(t, y, inputs)
         return states_eval - self.mass_matrix @ ydot
+
+
+class InitialConditions(SolverCallable):
+    "Returns initial conditions given inputs"
+
+    def __init__(self, function, model):
+        super().__init__(function, "initial conditions", model)
+        self.y_dummy = np.zeros(model.concatenated_initial_conditions.shape)
+
+    def __call__(self, inputs):
+        if self.form == "casadi":
+            if isinstance(inputs, dict):
+                inputs = casadi.vertcat(*[x for x in inputs.values()])
+            return self._function(0, self.y_dummy, inputs)
+        else:
+            return self._function(0, self.y_dummy, inputs=inputs).flatten()
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
index 70a5eda16..229cc7f6d 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -17,14 +17,18 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
     ----------
     tol : float, optional
         The tolerance for the solver (default is 1e-6).
+    extra_options : dict, optional
+        Any options to pass to the CasADi rootfinder.
+        Please consult `CasADi documentation <https://tinyurl.com/y7hrxm7d>`_ for
+        details.
     """
 
-    def __init__(self, method="lm", tol=1e-6, **extra_options):
+    def __init__(self, tol=1e-6, extra_options=None):
         super().__init__()
         self.tol = tol
         self.name = "CasADi algebraic solver"
         self.algebraic_solver = True
-        self.extra_options = extra_options
+        self.extra_options = extra_options or {}
         pybamm.citations.register("Andersson2019")
 
     @property
@@ -51,18 +55,32 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
             the solution will consist of `ProcessedSymbolicVariable` objects, which must
             be provided with inputs to obtain their value.
         """
-        y0 = model.y0
-
-        y = None
-
-        # Set up
         # Record whether there are any symbolic inputs
+        inputs = inputs or {}
         has_symbolic_inputs = any(isinstance(v, casadi.MX) for v in inputs.values())
 
         # Create casadi objects for the root-finder
         inputs = casadi.vertcat(*[x for x in inputs.values()])
+
+        y0 = model.y0
+        # The casadi algebraic solver can read rhs equations, but leaves them unchanged
+        # i.e. the part of the solution vector that corresponds to the differential
+        # equations will be equal to the initial condition provided. This allows this
+        # solver to be used for initialising the DAE solvers
+        if model.rhs == {}:
+            y0_diff = casadi.DM()
+            y0_alg = y0
+        else:
+            len_rhs = model.concatenated_rhs.size
+            y0_diff = y0[:len_rhs]
+            y0_alg = y0[len_rhs:]
+
+        y_alg = None
+
+        # Set up
         t_sym = casadi.MX.sym("t")
-        y_sym = casadi.MX.sym("y_alg", y0.shape[0])
+        y_alg_sym = casadi.MX.sym("y_alg", y0_alg.shape[0])
+        y_sym = casadi.vertcat(y0_diff, y_alg_sym)
         p_sym = casadi.MX.sym("p", inputs.shape[0])
 
         t_p_sym = casadi.vertcat(t_sym, p_sym)
@@ -72,7 +90,7 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
         roots = casadi.rootfinder(
             "roots",
             "newton",
-            dict(x=y_sym, p=t_p_sym, g=alg),
+            dict(x=y_alg_sym, p=t_p_sym, g=alg),
             {**self.extra_options, "abstol": self.tol},
         )
         for idx, t in enumerate(t_eval):
@@ -85,19 +103,20 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
                 pybamm.logger.debug(
                     "Keeping same solution at t={}".format(t * model.timescale_eval)
                 )
-                if y is None:
-                    y = y0
+                if y_alg is None:
+                    y_alg = y0_alg
                 else:
-                    y = casadi.horzcat(y, y0)
+                    y_alg = casadi.horzcat(y_alg, y0_alg)
             # Otherwise calculate new y_sol
             else:
                 t_inputs = casadi.vertcat(t, inputs)
                 # Solve
                 try:
-                    y_sol = roots(y0, t_inputs)
+                    y_alg_sol = roots(y0_alg, t_inputs)
                     success = True
                     message = None
                     # Check final output
+                    y_sol = casadi.vertcat(y0_diff, y_alg_sol)
                     fun = model.casadi_algebraic(t, y_sol, inputs)
                 except RuntimeError as err:
                     success = False
@@ -110,12 +129,12 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
                     has_symbolic_inputs is True or np.all(casadi.fabs(fun) < self.tol)
                 ):
                     # update initial guess for the next iteration
-                    y0 = y_sol
+                    y0_alg = y_alg_sol
                     # update solution array
-                    if y is None:
-                        y = y_sol
+                    if y_alg is None:
+                        y_alg = y_alg_sol
                     else:
-                        y = casadi.horzcat(y, y_sol)
+                        y_alg = casadi.horzcat(y_alg, y_alg_sol)
                 elif not success:
                     raise pybamm.SolverError(
                         "Could not find acceptable solution: {}".format(message)
@@ -131,5 +150,8 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
                         )
                     )
 
+        # Concatenate differential part
+        y_diff = casadi.horzcat(*[y0_diff] * len(t_eval))
+        y_sol = casadi.vertcat(y_diff, y_alg)
         # Return solution object (no events, so pass None to t_event, y_event)
-        return pybamm.Solution(t_eval, y, termination="success")
+        return pybamm.Solution(t_eval, y_sol, termination="success")
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index dd0bf0f35..cf45beb10 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -28,15 +28,24 @@ class CasadiSolver(pybamm.BaseSolver):
         The relative tolerance for the solver (default is 1e-6).
     atol : float, optional
         The absolute tolerance for the solver (default is 1e-6).
-    root_method : str, optional
-        The method to use for finding consistend initial conditions. Default is 'lm'.
+    root_method : str or pybamm algebraic solver class, optional
+        The method to use to find initial conditions (for DAE solvers).
+        If a solver class, must be an algebraic solver class.
+        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 root-finding. Default is 1e-6.
     max_step_decrease_counts : float, optional
         The maximum number of times step size can be decreased before an error is
         raised. Default is 5.
-    extra_options : keyword arguments, optional
-        Any extra keyword-arguments; these are passed directly to the CasADi integrator.
+    extra_options_setup : dict, optional
+        Any options to pass to the CasADi integrator when creating the integrator.
+        Please consult `CasADi documentation <https://tinyurl.com/y5rk76os>`_ for
+        details.
+    extra_options_call : dict, optional
+        Any options to pass to the CasADi integrator when calling the integrator.
         Please consult `CasADi documentation <https://tinyurl.com/y5rk76os>`_ for
         details.
 
@@ -50,7 +59,8 @@ class CasadiSolver(pybamm.BaseSolver):
         root_method="casadi",
         root_tol=1e-6,
         max_step_decrease_count=5,
-        **extra_options,
+        extra_options_setup=None,
+        extra_options_call=None,
     ):
         super().__init__("problem dependent", rtol, atol, root_method, root_tol)
         if mode in ["safe", "fast"]:
@@ -64,7 +74,8 @@ class CasadiSolver(pybamm.BaseSolver):
                 )
             )
         self.max_step_decrease_count = max_step_decrease_count
-        self.extra_options = extra_options
+        self.extra_options_setup = extra_options_setup or {}
+        self.extra_options_call = extra_options_call or {}
         self.name = "CasADi solver with '{}' mode".format(mode)
 
         # Initialize
@@ -91,13 +102,6 @@ class CasadiSolver(pybamm.BaseSolver):
         # convert inputs to casadi format
         inputs = casadi.vertcat(*[x for x in inputs.values()])
 
-        if len(model.rhs) == 0:
-            # casadi solver won't allow solving algebraic model so we have to raise an
-            # error here
-            raise pybamm.SolverError(
-                "Cannot use CasadiSolver to solve algebraic model, "
-                "use CasadiAlgebraicSolver instead"
-            )
         if self.mode == "fast":
             integrator = self.get_integrator(model, t_eval, inputs)
             solution = self._run_integrator(integrator, model, model.y0, inputs, t_eval)
@@ -110,18 +114,18 @@ class CasadiSolver(pybamm.BaseSolver):
             solution.termination = "final time"
             return solution
         elif self.mode == "safe":
+            y0 = model.y0
+            if isinstance(y0, casadi.DM):
+                y0 = y0.full().flatten()
             # Step-and-check
             t = t_eval[0]
             init_event_signs = np.sign(
                 np.concatenate(
-                    [
-                        event(t, model.y0, inputs)
-                        for event in model.terminate_events_eval
-                    ]
+                    [event(t, y0, inputs) for event in model.terminate_events_eval]
                 )
             )
             pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))
-            y0 = model.y0
+
             # Initialize solution
             solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis])
             solution.solve_time = 0
@@ -189,11 +193,11 @@ class CasadiSolver(pybamm.BaseSolver):
             algebraic = model.casadi_algebraic
 
             options = {
+                **self.extra_options_setup,
                 "grid": t_eval,
                 "reltol": self.rtol,
                 "abstol": self.atol,
                 "output_t0": True,
-                "max_num_steps": self.max_steps,
             }
 
             # set up and solve
@@ -228,11 +232,11 @@ class CasadiSolver(pybamm.BaseSolver):
         )
 
     def _run_integrator(self, integrator, model, y0, inputs, t_eval):
-        rhs_size = model.rhs_eval(t_eval[0], y0, inputs).shape[0]
+        rhs_size = model.concatenated_rhs.size
         y0_diff, y0_alg = np.split(y0, [rhs_size])
         try:
             # Try solving
-            sol = integrator(x0=y0_diff, z0=y0_alg, p=inputs, **self.extra_options)
+            sol = integrator(x0=y0_diff, z0=y0_alg, p=inputs, **self.extra_options_call)
             y_values = np.concatenate([sol["xf"].full(), sol["zf"].full()])
             return pybamm.Solution(t_eval, y_values)
         except RuntimeError as e:
diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py
index 8909a78bb..d90a519ed 100644
--- a/pybamm/solvers/idaklu_solver.py
+++ b/pybamm/solvers/idaklu_solver.py
@@ -27,17 +27,24 @@ class IDAKLUSolver(pybamm.BaseSolver):
         The relative tolerance for the solver (default is 1e-6).
     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")
+    root_method : str or pybamm algebraic solver class, optional
+        The method to use to find initial conditions (for DAE solvers).
+        If a solver class, must be an algebraic solver class.
+        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-8).
-    max_steps: int, optional
-        The maximum number of steps the solver will take before terminating
-        (default is 1000).
     """
 
     def __init__(
-        self, rtol=1e-6, atol=1e-6, root_method="casadi", root_tol=1e-6, max_steps=1000
+        self,
+        rtol=1e-6,
+        atol=1e-6,
+        root_method="casadi",
+        root_tol=1e-6,
+        max_steps="deprecated",
     ):
 
         if idaklu_spec is None:
@@ -158,9 +165,13 @@ class IDAKLUSolver(pybamm.BaseSolver):
         except AttributeError:
             atol = self._atol
 
-        rtol = self._rtol
-        atol = self._check_atol_type(atol, model.y0.size)
         y0 = model.y0
+        if isinstance(y0, casadi.DM):
+            y0 = y0.full().flatten()
+
+        rtol = self._rtol
+        atol = self._check_atol_type(atol, y0.size)
+
         mass_matrix = model.mass_matrix.entries
 
         if model.jacobian_eval:
@@ -242,7 +253,7 @@ class IDAKLUSolver(pybamm.BaseSolver):
 
         t = sol.t
         number_of_timesteps = t.size
-        number_of_states = model.y0.size
+        number_of_states = y0.size
         y_out = sol.y.reshape((number_of_timesteps, number_of_states))
 
         # return solution, we need to tranpose y to match scipy's interface
diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py
index 77adbe3dd..6f3b56ec8 100644
--- a/pybamm/solvers/scikits_dae_solver.py
+++ b/pybamm/solvers/scikits_dae_solver.py
@@ -27,13 +27,22 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
         The relative tolerance for the solver (default is 1e-6).
     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")
+    root_method : str or pybamm algebraic solver class, optional
+        The method to use to find initial conditions (for DAE solvers).
+        If a solver class, must be an algebraic solver class.
+        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
-        The maximum number of steps the solver will take before terminating
-        (default is 1000).
+    extra_options : dict, optional
+        Any options to pass to the solver.
+        Please consult `scikits.odes documentation
+        <https://bmcage.github.io/odes/dev/index.html>`_ for details.
+        Some common keys:
+
+        - 'max_steps': maximum (int) number of steps the solver can take
     """
 
     def __init__(
@@ -43,7 +52,8 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
         atol=1e-6,
         root_method="casadi",
         root_tol=1e-6,
-        max_steps=1000,
+        extra_options=None,
+        max_steps="deprecated",
     ):
         if scikits_odes_spec is None:
             raise ImportError("scikits.odes is not installed")
@@ -51,6 +61,8 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
         super().__init__(method, rtol, atol, root_method, root_tol, max_steps)
         self.name = "Scikits DAE solver ({})".format(method)
 
+        self.extra_options = extra_options or {}
+
         pybamm.citations.register("scikits-odes")
         pybamm.citations.register("hindmarsh2000pvode")
         pybamm.citations.register("hindmarsh2005sundials")
@@ -72,8 +84,11 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
         if model.convert_to_format == "casadi":
             inputs = casadi.vertcat(*[x for x in inputs.values()])
 
-        residuals = model.residuals_eval
         y0 = model.y0
+        if isinstance(y0, casadi.DM):
+            y0 = y0.full().flatten()
+
+        residuals = model.residuals_eval
         events = model.terminate_events_eval
         jacobian = model.jacobian_eval
         mass_matrix = model.mass_matrix.entries
@@ -85,10 +100,10 @@ class ScikitsDaeSolver(pybamm.BaseSolver):
             return_root[:] = [event(t, y, inputs) for event in events]
 
         extra_options = {
+            **self.extra_options,
             "old_api": False,
             "rtol": self.rtol,
             "atol": self.atol,
-            "max_steps": self.max_steps,
         }
 
         if jacobian:
diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py
index 21536420f..457e5b520 100644
--- a/pybamm/solvers/scikits_ode_solver.py
+++ b/pybamm/solvers/scikits_ode_solver.py
@@ -31,16 +31,34 @@ class ScikitsOdeSolver(pybamm.BaseSolver):
         The relative tolerance for the solver (default is 1e-6).
     atol : float, optional
         The absolute tolerance for the solver (default is 1e-6).
-    linsolver : str, optional
-            Can be 'dense' (= default), 'lapackdense', 'spgmr', 'spbcgs', 'sptfqmr'
+    extra_options : dict, optional
+        Any options to pass to the solver.
+        Please consult `scikits.odes documentation
+        <https://bmcage.github.io/odes/dev/index.html>`_ for details.
+        Some common keys:
+
+        - 'linsolver': can be 'dense' (= default), 'lapackdense', 'spgmr', 'spbcgs', \
+        'sptfqmr'
     """
 
-    def __init__(self, method="cvode", rtol=1e-6, atol=1e-6, linsolver="dense"):
+    def __init__(
+        self,
+        method="cvode",
+        rtol=1e-6,
+        atol=1e-6,
+        linsolver="deprecated",
+        extra_options=None,
+    ):
         if scikits_odes_spec is None:
             raise ImportError("scikits.odes is not installed")
 
         super().__init__(method, rtol, atol)
-        self.linsolver = linsolver
+        self.extra_options = extra_options or {}
+        if linsolver != "deprecated":
+            raise ValueError(
+                "linsolver has been deprecated. Pass 'linsolver' to extra_options "
+                "dictionary instead"
+            )
         self.ode_solver = True
         self.name = "Scikits ODE solver ({})".format(method)
 
@@ -65,8 +83,11 @@ class ScikitsOdeSolver(pybamm.BaseSolver):
         if model.rhs_eval.form == "casadi":
             inputs = casadi.vertcat(*[x for x in inputs.values()])
 
-        derivs = model.rhs_eval
         y0 = model.y0
+        if isinstance(y0, casadi.DM):
+            y0 = y0.full().flatten()
+
+        derivs = model.rhs_eval
         events = model.terminate_events_eval
         jacobian = model.jacobian_eval
 
@@ -101,16 +122,19 @@ class ScikitsOdeSolver(pybamm.BaseSolver):
                 return 0
 
         extra_options = {
+            **self.extra_options,
             "old_api": False,
             "rtol": self.rtol,
             "atol": self.atol,
-            "linsolver": self.linsolver,
         }
 
+        # Read linsolver (defaults to dense)
+        linsolver = extra_options.get("linsolver", "dense")
+
         if jacobian:
-            if self.linsolver in ("dense", "lapackdense"):
+            if linsolver in ("dense", "lapackdense"):
                 extra_options.update({"jacfn": jacfn})
-            elif self.linsolver in ("spgmr", "spbcgs", "sptfqmr"):
+            elif linsolver in ("spgmr", "spbcgs", "sptfqmr"):
                 extra_options.update(
                     {
                         "jac_times_setupfn": jac_times_setupfn,
diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py
index 45389f146..613ae72a5 100644
--- a/pybamm/solvers/scipy_solver.py
+++ b/pybamm/solvers/scipy_solver.py
@@ -9,7 +9,7 @@ import numpy as np
 
 
 class ScipySolver(pybamm.BaseSolver):
-    """Solve a discretised model, using scipy._integrate.solve_ivp.
+    """Solve a discretised model, using scipy.integrate.solve_ivp.
 
     Parameters
     ----------
@@ -19,11 +19,16 @@ class ScipySolver(pybamm.BaseSolver):
         The relative tolerance for the solver (default is 1e-6).
     atol : float, optional
         The absolute tolerance for the solver (default is 1e-6).
+    extra_options : dict, optional
+        Any options to pass to the solver.
+        Please consult `SciPy documentation <https://tinyurl.com/yafgqg9y>`_ for
+        details.
     """
 
-    def __init__(self, method="BDF", rtol=1e-6, atol=1e-6):
+    def __init__(self, method="BDF", rtol=1e-6, atol=1e-6, extra_options=None):
         super().__init__(method, rtol, atol)
         self.ode_solver = True
+        self.extra_options = extra_options or {}
         self.name = "Scipy solver ({})".format(method)
         pybamm.citations.register("virtanen2020scipy")
 
@@ -50,7 +55,12 @@ class ScipySolver(pybamm.BaseSolver):
         if model.convert_to_format == "casadi":
             inputs = casadi.vertcat(*[x for x in inputs.values()])
 
-        extra_options = {"rtol": self.rtol, "atol": self.atol}
+        extra_options = {**self.extra_options, "rtol": self.rtol, "atol": self.atol}
+
+        # Initial conditions
+        y0 = model.y0
+        if isinstance(y0, casadi.DM):
+            y0 = y0.full().flatten()
 
         # check for user-supplied Jacobian
         implicit_methods = ["Radau", "BDF", "LSODA"]
@@ -76,7 +86,7 @@ class ScipySolver(pybamm.BaseSolver):
         sol = it.solve_ivp(
             lambda t, y: model.rhs_eval(t, y, inputs),
             (t_eval[0], t_eval[-1]),
-            model.y0,
+            y0,
             t_eval=t_eval,
             method=self.method,
             dense_output=True,

Pass extra options to solvers
it would be good to optionally get some extra solver output, e.g. from casadi we can readily get something like 
```
Number of steps taken by SUNDIALS: 394
Number of calls to the user’s f function: 716
Number of calls made to the linear solver setup function: 39
Number of error test failures: 4
Method order used on the last internal step: 3
Method order to be used on the next internal step: 3
Actual value of initial step size: 2.56788e-08
Step size taken on the last internal step: 0.000225891
Step size to be attempted on the next internal step: 0.000225891
Current internal time reached: 0.000225891
Number of nonlinear iterations performed: 711
Number of nonlinear convergence failures: 3
```


diff --git a/tests/unit/test_solvers/test_algebraic_solver.py b/tests/unit/test_solvers/test_algebraic_solver.py
index b4ef18ecf..45c858327 100644
--- a/tests/unit/test_solvers/test_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_algebraic_solver.py
@@ -9,8 +9,11 @@ from tests import get_discretisation_for_testing
 
 class TestAlgebraicSolver(unittest.TestCase):
     def test_algebraic_solver_init(self):
-        solver = pybamm.AlgebraicSolver(method="hybr", tol=1e-4)
+        solver = pybamm.AlgebraicSolver(
+            method="hybr", tol=1e-4, extra_options={"maxfev": 100}
+        )
         self.assertEqual(solver.method, "hybr")
+        self.assertEqual(solver.extra_options, {"maxfev": 100})
         self.assertEqual(solver.tol, 1e-4)
 
         solver.method = "krylov"
@@ -37,21 +40,31 @@ class TestAlgebraicSolver(unittest.TestCase):
         # Simple system: a single algebraic equation
         class Model:
             y0 = np.array([2])
-            jacobian_eval = None
+            rhs = {}
+            timescale_eval = 1
+            jac_algebraic_eval = None
             convert_to_format = "python"
 
             def algebraic_eval(self, t, y, inputs):
                 return y + 2
 
-        solver = pybamm.AlgebraicSolver()
+        # Try passing extra options to solver
+        solver = pybamm.AlgebraicSolver(extra_options={"maxiter": 100})
         model = Model()
         solution = solver._integrate(model, np.array([0]))
         np.testing.assert_array_equal(solution.y, -2)
 
+        # Relax options and see worse results
+        solver = pybamm.AlgebraicSolver(extra_options={"ftol": 1})
+        solution = solver._integrate(model, np.array([0]))
+        self.assertNotEqual(solution.y, -2)
+
     def test_root_find_fail(self):
         class Model:
             y0 = np.array([2])
-            jacobian_eval = None
+            rhs = {}
+            timescale_eval = 1
+            jac_algebraic_eval = None
             convert_to_format = "casadi"
 
             def algebraic_eval(self, t, y, inputs):
@@ -79,12 +92,14 @@ class TestAlgebraicSolver(unittest.TestCase):
 
         class Model:
             y0 = np.zeros(2)
+            rhs = {}
+            timescale_eval = 1
             convert_to_format = "python"
 
             def algebraic_eval(self, t, y, inputs):
                 return A @ y - b
 
-            def jacobian_eval(self, t, y, inputs):
+            def jac_algebraic_eval(self, t, y, inputs):
                 return A
 
         model = Model()
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 0540a874d..885032e5c 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -20,6 +20,27 @@ class TestBaseSolver(unittest.TestCase):
         solver.rtol = 1e-7
         self.assertEqual(solver.rtol, 1e-7)
 
+        # max_steps deprecated
+        with self.assertRaisesRegex(ValueError, "max_steps has been deprecated"):
+            pybamm.BaseSolver(max_steps=10)
+
+    def test_root_method_init(self):
+        solver = pybamm.BaseSolver(root_method="casadi")
+        self.assertIsInstance(solver.root_method, pybamm.CasadiAlgebraicSolver)
+
+        solver = pybamm.BaseSolver(root_method="lm")
+        self.assertIsInstance(solver.root_method, pybamm.AlgebraicSolver)
+        self.assertEqual(solver.root_method.method, "lm")
+
+        root_solver = pybamm.AlgebraicSolver()
+        solver = pybamm.BaseSolver(root_method=root_solver)
+        self.assertEqual(solver.root_method, root_solver)
+
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Root method must be an algebraic solver"
+        ):
+            pybamm.BaseSolver(root_method=pybamm.ScipySolver())
+
     def test_step_or_solve_empty_model(self):
         model = pybamm.BaseModel()
         solver = pybamm.BaseSolver()
@@ -74,9 +95,10 @@ class TestBaseSolver(unittest.TestCase):
         # Simple system: a single algebraic equation
         class ScalarModel:
             def __init__(self):
-                self.concatenated_initial_conditions = np.array([[2]])
+                self.y0 = np.array([2])
+                self.rhs = {}
                 self.jac_algebraic_eval = None
-                self.timescale = 1
+                self.timescale_eval = 1
                 t = casadi.MX.sym("t")
                 y = casadi.MX.sym("y")
                 p = casadi.MX.sym("p")
@@ -106,9 +128,11 @@ class TestBaseSolver(unittest.TestCase):
 
         class VectorModel:
             def __init__(self):
-                self.concatenated_initial_conditions = np.zeros_like(vec)
+                self.y0 = np.zeros_like(vec)
+                self.rhs = {"test": "test"}
+                self.concatenated_rhs = np.array([1])
                 self.jac_algebraic_eval = None
-                self.timescale = 1
+                self.timescale_eval = 1
                 t = casadi.MX.sym("t")
                 y = casadi.MX.sym("y", vec.size)
                 p = casadi.MX.sym("p")
@@ -151,9 +175,10 @@ class TestBaseSolver(unittest.TestCase):
     def test_fail_consistent_initial_conditions(self):
         class Model:
             def __init__(self):
-                self.concatenated_initial_conditions = np.array([2])
+                self.y0 = np.array([2])
+                self.rhs = {}
                 self.jac_algebraic_eval = None
-                self.timescale = 1
+                self.timescale_eval = 1
                 t = casadi.MX.sym("t")
                 y = casadi.MX.sym("y")
                 p = casadi.MX.sym("p")
@@ -173,20 +198,18 @@ class TestBaseSolver(unittest.TestCase):
 
         with self.assertRaisesRegex(
             pybamm.SolverError,
-            "Could not find consistent initial conditions: The iteration is not making",
+            "Could not find acceptable solution: The iteration is not making",
         ):
             solver.calculate_consistent_state(Model())
         solver = pybamm.BaseSolver(root_method="lm")
         with self.assertRaisesRegex(
-            pybamm.SolverError,
-            "Could not find consistent initial conditions: solver terminated",
+            pybamm.SolverError, "Could not find acceptable solution: 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",
+            pybamm.SolverError, "Could not find acceptable solution: .../casadi",
         ):
             solver.calculate_consistent_state(Model())
 
@@ -224,7 +247,7 @@ class TestBaseSolver(unittest.TestCase):
         disc = pybamm.Discretisation()
         disc.process_model(model)
 
-        solver = pybamm.BaseSolver()
+        solver = pybamm.BaseSolver(root_method="casadi")
         pybamm.set_logging_level("ERROR")
         solver.set_up(model, {})
         self.assertEqual(model.convert_to_format, "casadi")
diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index efe1c0bab..d9fc2531e 100644
--- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -56,6 +56,7 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
             t = casadi.MX.sym("t")
             y = casadi.MX.sym("y")
             p = casadi.MX.sym("p")
+            rhs = {}
             casadi_algebraic = casadi.Function("alg", [t, y, p], [y ** 2 + 1])
 
             def algebraic_eval(self, t, y, inputs):
@@ -69,7 +70,7 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
             pybamm.SolverError, "Could not find acceptable solution: .../casadi"
         ):
             solver._integrate(model, np.array([0]), {})
-        solver = pybamm.CasadiAlgebraicSolver(error_on_fail=False)
+        solver = pybamm.CasadiAlgebraicSolver(extra_options={"error_on_fail": False})
         with self.assertRaisesRegex(
             pybamm.SolverError, "Could not find acceptable solution: solver terminated"
         ):
@@ -225,6 +226,51 @@ class TestCasadiAlgebraicSolverSensitivity(unittest.TestCase):
             solution["var"].sensitivity({"param": p}), -np.eye(40)
         )
 
+    def test_solve_with_symbolic_input_in_initial_conditions(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var")
+        model = pybamm.BaseModel()
+        model.algebraic = {var: var + 2}
+        model.initial_conditions = {var: pybamm.InputParameter("param")}
+        model.variables = {"var": var}
+
+        # create discretisation
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, [0])
+        np.testing.assert_array_equal(solution["var"].value({"param": 7}), -2)
+        np.testing.assert_array_equal(solution["var"].value({"param": 3}), -2)
+        np.testing.assert_array_equal(solution["var"].sensitivity({"param": 3}), 0)
+
+    def test_least_squares_fit_input_in_initial_conditions(self):
+        # Simple system: a single algebraic equation
+        var = pybamm.Variable("var", domain="negative electrode")
+        model = pybamm.BaseModel()
+        p = pybamm.InputParameter("p")
+        q = pybamm.InputParameter("q")
+        model.algebraic = {var: (var - p)}
+        model.initial_conditions = {var: p}
+        model.variables = {"objective": (var - q) ** 2 + (p - 3) ** 2}
+
+        # create discretisation
+        disc = tests.get_discretisation_for_testing()
+        disc.process_model(model)
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver()
+        solution = solver.solve(model, [0])
+        sol_var = solution["objective"]
+
+        def objective(x):
+            return sol_var.value({"p": x[0], "q": x[1]}).full().flatten()
+
+        # without jacobian
+        lsq_sol = least_squares(objective, [2, 2], method="lm")
+        np.testing.assert_array_almost_equal(lsq_sol.x, [3, 3], decimal=3)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 1018baa68..b9ee3b613 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -81,7 +81,7 @@ class TestCasadiSolver(unittest.TestCase):
         disc = pybamm.Discretisation()
         disc.process_model(model)
 
-        solver = pybamm.CasadiSolver(regularity_check=False)
+        solver = pybamm.CasadiSolver(extra_options_call={"regularity_check": False})
 
         # Solve with failure at t=2
         t_eval = np.linspace(0, 20, 100)
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index c646075a3..1bab2b4d4 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -11,6 +11,11 @@ import sys
 
 @unittest.skipIf(not pybamm.have_scikits_odes(), "scikits.odes not installed")
 class TestScikitsSolvers(unittest.TestCase):
+    def test_init(self):
+        # linsolver deprecated
+        with self.assertRaisesRegex(ValueError, "linsolver has been deprecated"):
+            pybamm.ScikitsOdeSolver(linsolver="lapackdense")
+
     def test_model_ode_integrate_failure(self):
         # Turn off warnings to ignore sqrt error
         warnings.simplefilter("ignore")
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index e7393927c..fd21be080 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -26,7 +26,10 @@ class TestScipySolver(unittest.TestCase):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
         # Solve
-        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8, method="RK45")
+        # Make sure that passing in extra options works
+        solver = pybamm.ScipySolver(
+            rtol=1e-8, atol=1e-8, method="RK45", extra_options={"first_step": 1e-4}
+        )
         t_eval = np.linspace(0, 1, 80)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gcc gfortran libopenblas-dev liblapack-dev wget
3.9
requirements.txt requirements-docs.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.16
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
docutils==0.21.2
exceptiongroup==1.2.2
executing==2.2.0
fastjsonschema==2.21.1
flake8==7.2.0
fonttools==4.56.0
fqdn==1.5.1
guzzle_sphinx_theme==0.7.11
h11==0.14.0
httpcore==1.0.7
httpx==0.28.1
idna==3.10
imagesize==1.4.1
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
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@23826728806b50c10b8197268a5e513289934392#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
pytz==2025.2
PyYAML==6.0.2
pyzmq==26.3.0
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
snowballstemmer==2.2.0
soupsieve==2.6
Sphinx==7.4.7
sphinxcontrib-applehelp==2.0.0
sphinxcontrib-devhelp==2.0.0
sphinxcontrib-htmlhelp==2.1.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==2.0.0
sphinxcontrib-serializinghtml==2.0.0
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:
      - alabaster==0.7.16
      - 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
      - docutils==0.21.2
      - exceptiongroup==1.2.2
      - executing==2.2.0
      - fastjsonschema==2.21.1
      - flake8==7.2.0
      - fonttools==4.56.0
      - fqdn==1.5.1
      - guzzle-sphinx-theme==0.7.11
      - h11==0.14.0
      - httpcore==1.0.7
      - httpx==0.28.1
      - idna==3.10
      - imagesize==1.4.1
      - 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
      - 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
      - pytz==2025.2
      - pyyaml==6.0.2
      - pyzmq==26.3.0
      - 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
      - snowballstemmer==2.2.0
      - soupsieve==2.6
      - sphinx==7.4.7
      - sphinxcontrib-applehelp==2.0.0
      - sphinxcontrib-devhelp==2.0.0
      - sphinxcontrib-htmlhelp==2.1.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==2.0.0
      - sphinxcontrib-serializinghtml==2.0.0
      - 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_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_algebraic_solver_init tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_root_find_fail tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_simple_root_find tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_with_jacobian tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_base_solver_init tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_find_consistent_initial_conditions tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_root_method_init tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_least_squares_fit_input_in_initial_conditions tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input_in_initial_conditions tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_python
tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_fail_consistent_initial_conditions tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_root_find_fail
tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver_with_time tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_solve_with_input tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_wrong_solver tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_block_symbolic_inputs tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_convert_to_casadi_format tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_discretise_model 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_t_eval_none tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_algebraic_solver_init tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_model_solver_with_time tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_simple_root_find tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_simple_root_find_correct_initial_guess tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_solve_with_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_least_squares_fit tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input_1D_scalar_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input_1D_vector_input tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_bad_mode tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_dae_solver_algebraic_model tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_dae_inputs_in_initial_conditions tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_events tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_failure tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_python tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_with_external tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_with_inputs tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_solver_with_non_identity_mass tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_step tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_step_events tests/unit/test_solvers/test_casadi_solver.py::TestCasadiSolver::test_model_step_with_input 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_inputs_in_initial_conditions 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_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/scikits_ode_solver.py pybamm/solvers/casadi_solver.py pybamm/solvers/scipy_solver.py pybamm/solvers/casadi_algebraic_solver.py pybamm/solvers/algebraic_solver.py pybamm/solvers/idaklu_solver.py pybamm/solvers/base_solver.py pybamm/solvers/scikits_dae_solver.py
pybamm/solvers/scikits_ode_solver.py pybamm/solvers/casadi_solver.py pybamm/solvers/scipy_solver.py pybamm/solvers/casadi_algebraic_solver.py pybamm/solvers/algebraic_solver.py pybamm/solvers/idaklu_solver.py pybamm/solvers/base_solver.py pybamm/solvers/scikits_dae_solver.py