function varargout = ode_bdf(odefun, tspan, y0, options)

cvode_options = CVodeSetOptions(...
    'LMM', 'BDF',...
    'NonlinearSolver', 'Newton',...
    'RelTol', options.RelTol,...
);