def integrate(self, x, y0, params=(), atol=1e-8, rtol=1e-8, **kwargs):
        """ Integrate the system of ordinary differential equations.

        Solves the initial value problem (IVP).

        Parameters
        ----------
        x : array_like or pair (start and final time) or float
            if float:
                make it a pair: (0, x)
            if pair or length-2 array:
                initial and final value of the independent variable
            if array_like:
                values of independent variable report at
        y0 : array_like
            Initial values at x[0] for the dependent variables.
        params : array_like (default: tuple())
            Value of parameters passed to user-supplied callbacks.
        integrator : str or None
            Name of integrator, one of:
                - 'scipy': :meth:`_integrate_scipy`
                - 'gsl': :meth:`_integrate_gsl`
                - 'odeint': :meth:`_integrate_odeint`
                - 'cvode':  :meth:`_integrate_cvode`

            See respective method for more information.
            If ``None``: ``os.environ.get('PYODESYS_INTEGRATOR', 'scipy')``
        atol : float
            Absolute tolerance
        rtol : float
            Relative tolerance
        with_jacobian : bool or None (default)
            Whether to use the jacobian. When ``None`` the choice is
            done automatically (only used when required). This matters
            when jacobian is derived at runtime (high computational cost).
        with_jtimes : bool (default: False)
            Whether to use the jacobian-vector product. This is only supported
            by ``cvode`` and only when ``linear_solver`` is one of: gmres',
            'gmres_classic', 'bicgstab', 'tfqmr'. See the documentation
            for ``pycvodes`` for more information.
        force_predefined : bool (default: False)
            override behaviour of ``len(x) == 2`` => :meth:`adaptive`
        \\*\\*kwargs :
            Additional keyword arguments for ``_integrate_$(integrator)``.

        Returns
        -------
        Length 3 tuple: (x, yout, info)
            x : array of values of the independent variable
            yout : array of the dependent variable(s) for the different
                values of x.
            info : dict ('nfev' is guaranteed to be a key)
        """
        arrs = self.to_arrays(x, y0, params)
        _x, _y, _p = _arrs = self.pre_process(*arrs)
        ndims = [a.ndim for a in _arrs]
        if ndims == [1, 1, 1]:
            twodim = False
        elif ndims == [2, 2, 2]:
            twodim = True
        else:
            raise ValueError("Pre-processor made ndims inconsistent?")

        if self.append_iv:
            _p = self.numpy.concatenate((_p, _y), axis=-1)

        if hasattr(self, 'ny'):
            if _y.shape[-1] != self.ny:
                raise ValueError("Incorrect shape of intern_y0")
        if isinstance(atol, dict):
            kwargs['atol'] = [atol[k] for k in self.names]
        else:
            kwargs['atol'] = atol
        kwargs['rtol'] = rtol

        integrator = kwargs.pop('integrator', None)
        if integrator is None:
            integrator = os.environ.get('PYODESYS_INTEGRATOR', 'scipy')

        args = tuple(map(self.numpy.atleast_2d, (_x, _y, _p)))

        self._current_integration_kwargs = kwargs
        if isinstance(integrator, str):
            nfo = getattr(self, '_integrate_' + integrator)(*args, **kwargs)
        else:
            kwargs['with_jacobian'] = getattr(integrator, 'with_jacobian', None)
            nfo = self._integrate(integrator.integrate_adaptive,
                                  integrator.integrate_predefined,
                                  *args, **kwargs)
        if twodim:
            _xout = [d['internal_xout'] for d in nfo]
            _yout = [d['internal_yout'] for d in nfo]
            _params = [d['internal_params'] for d in nfo]
            res = [Result(*(self.post_process(_xout[i], _yout[i], _params[i]) + (nfo[i], self)))
                   for i in range(len(nfo))]
        else:
            _xout = nfo[0]['internal_xout']
            _yout = nfo[0]['internal_yout']

            self._internal = _xout.copy(), _yout.copy(), _p.copy()
            nfo = nfo[0]
            res = Result(*(self.post_process(_xout, _yout, _p) + (nfo, self)))
        return res
Integrate the system of ordinary differential equations.

        Solves the initial value problem (IVP).

        Parameters
        ----------
        x : array_like or pair (start and final time) or float
            if float:
                make it a pair: (0, x)
            if pair or length-2 array:
                initial and final value of the independent variable
            if array_like:
                values of independent variable report at
        y0 : array_like
            Initial values at x[0] for the dependent variables.
        params : array_like (default: tuple())
            Value of parameters passed to user-supplied callbacks.
        integrator : str or None
            Name of integrator, one of:
                - 'scipy': :meth:`_integrate_scipy`
                - 'gsl': :meth:`_integrate_gsl`
                - 'odeint': :meth:`_integrate_odeint`
                - 'cvode':  :meth:`_integrate_cvode`

            See respective method for more information.
            If ``None``: ``os.environ.get('PYODESYS_INTEGRATOR', 'scipy')``
        atol : float
            Absolute tolerance
        rtol : float
            Relative tolerance
        with_jacobian : bool or None (default)
            Whether to use the jacobian. When ``None`` the choice is
            done automatically (only used when required). This matters
            when jacobian is derived at runtime (high computational cost).
        with_jtimes : bool (default: False)
            Whether to use the jacobian-vector product. This is only supported
            by ``cvode`` and only when ``linear_solver`` is one of: gmres',
            'gmres_classic', 'bicgstab', 'tfqmr'. See the documentation
            for ``pycvodes`` for more information.
        force_predefined : bool (default: False)
            override behaviour of ``len(x) == 2`` => :meth:`adaptive`
        \\*\\*kwargs :
            Additional keyword arguments for ``_integrate_$(integrator)``.

        Returns
        -------
        Length 3 tuple: (x, yout, info)
            x : array of values of the independent variable
            yout : array of the dependent variable(s) for the different
                values of x.
            info : dict ('nfev' is guaranteed to be a key)