def METHOD_NAME(self, model, t_eval, inputs_dict=None):
    """
    Solve a model defined by dydt with initial conditions y0.
    Parameters
    ----------
    model : :class:`pybamm.BaseModel`
        The model whose solution to calculate.
    t_eval : numeric type
        The times at which to compute the solution
    inputs_dict : dict, optional
        Any input parameters to pass to the model when solving
    """
    inputs_dict = inputs_dict or {}
    if model.convert_to_format == "casadi":
        inputs = casadi.vertcat(*[x for x in inputs_dict.values()])
    else:
        inputs = inputs_dict
    y0 = model.y0
    if isinstance(y0, casadi.DM):
        y0 = y0.full()
    y0 = y0.flatten()
    rhs_algebraic_eval = model.rhs_algebraic_eval
    events = model.terminate_events_eval
    jacobian = model.jac_rhs_algebraic_eval
    if model.convert_to_format == "jax":
        mass_matrix = model.mass_matrix.entries.toarray()
    else:
        mass_matrix = model.mass_matrix.entries
    if model.convert_to_format == "casadi":
        def eqsres(t, y, ydot, return_residuals):
            return_residuals[:] = (
                rhs_algebraic_eval(t, y, inputs).full().flatten()
                - mass_matrix @ ydot
            )
    else:
        def eqsres(t, y, ydot, return_residuals):
            return_residuals[:] = (
                rhs_algebraic_eval(t, y, inputs).flatten() - mass_matrix @ ydot
            )
    def rootfn(t, y, ydot, return_root):
        return_root[:] = [float(event(t, y, inputs)) for event in events]
    extra_options = {
        **self.extra_options,
        "old_api": False,
        "rtol": self.rtol,
        "atol": self.atol,
    }
    if jacobian:
        jac_y0_t0 = jacobian(t_eval[0], y0, inputs)
        if sparse.issparse(jac_y0_t0):
            def jacfn(t, y, ydot, residuals, cj, J):
                jac_eval = jacobian(t, y, inputs) - cj * mass_matrix
                J[:][:] = jac_eval.toarray()
        else:
            def jacfn(t, y, ydot, residuals, cj, J):
                jac_eval = jacobian(t, y, inputs) - cj * mass_matrix
                J[:][:] = jac_eval
        extra_options.update({"jacfn": jacfn})
    if events:
        extra_options.update({"rootfn": rootfn, "nr_rootfns": len(events)})
    # solver works with ydot0 set to zero
    ydot0 = np.zeros_like(y0)
    # set up and solve
    dae_solver = scikits_odes.dae(self.method, eqsres, **extra_options)
    timer = pybamm.Timer()
    sol = dae_solver.solve(t_eval, y0, ydot0)
    integration_time = timer.time()
    # return solution, we need to tranpose y to match scipy's interface
    if sol.flag in [0, 2]:
        # 0 = solved for all t_eval
        if sol.flag == 0:
            termination = "final time"
        # 2 = found root(s)
        elif sol.flag == 2:
            termination = "event"
        if sol.roots.t is None:
            t_root = None
        else:
            t_root = sol.roots.t
        sol = pybamm.Solution(
            sol.values.t,
            np.transpose(sol.values.y),
            model,
            inputs_dict,
            t_root,
            np.transpose(sol.roots.y),
            termination,
        )
        sol.integration_time = integration_time
        return sol
    else:
        raise pybamm.SolverError(sol.message)