code: def integrate(self, x, y0, params=(), atol=1e-8, rtol=1e-8, **kwargs):
          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