How to use the pyodesys.results.Result function in pyodesys

def plot_reaction_contributions(
        xyp, rsys, rate_exprs_cb, substance_keys=None, varied=None, axes=None,
        total=False, linthreshy=1e-9, relative=False, xscale='log', yscale='symlog',
        xlabel='Time', ylabel=None, combine_equilibria=False, selection=slice(None),
        unit_registry=None):
    """ Plots per reaction contributions to concentration evolution of a substance.

    Parameters
    ----------
    xyp : ``pyodesys.results.Result`` instance or length 3 tuple or xout,yout,params
    result : pyodesys.results.Result
    substance_key : str

    """
    from pyodesys.results import Result
    if isinstance(xyp, Result):
        xyp = xyp.odesys.to_arrays(xyp.xout, xyp.yout, xyp.params, reshape=False)
    if varied is None:
        varied = xyp[0]
    if xyp[1].shape[-2] != varied.size:
        raise ValueError("Size mismatch between varied and yout")
    if substance_keys is None:
        substance_keys = rsys.substances.keys()
    if axes is None:
        _fig, axes = plt.subplots(len(substance_keys))
    rates = rate_exprs_cb(*xyp)
    if unit_registry is not None:
        time_unit = get_derived_unit(unit_registry, 'time')
        conc_unit = get_derived_unit(unit_registry, 'concentration')
        rates = to_unitless(rates*conc_unit/time_unit, u.molar/u.second)

    eqk1, eqk2, eqs = _combine_rxns_to_eq(rsys) if combine_equilibria else ([], [], [])

def copy(self):
        return Result(self.xout.copy(), self.yout.copy(), self.params.copy(),
                      self.info.copy(), self.odesys)

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

if multimode:
                _x, y0 = [], []
                for idx in range(multimode):
                    _x.append(_new_x(res[idx].xout, x[idx], next_autonomous))
                    y0.append(res[idx].yout[-1, :])
                    if next_autonomous:
                        glob_x[idx] += res[idx].xout[-1]
                x = _x
            else:
                x = _new_x(res.xout, x, next_autonomous)
                y0 = res.yout[-1, :]
                if next_autonomous:
                    glob_x += res.xout[-1]
    if multimode:  # don't return defaultdict
        tot_nfo = [dict(nsys=oi+1, **_nfo) for _nfo in tot_nfo]
        return [Result(tot_x[idx], tot_y[idx], res[idx].params, tot_nfo[idx], odes[0])
                for idx in range(len(res))]
    else:
        tot_nfo = dict(nsys=oi+1, **tot_nfo)
        return Result(tot_x, tot_y, res.params, tot_nfo, odes[0])

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

y0.append(res[idx].yout[-1, :])
                    if next_autonomous:
                        glob_x[idx] += res[idx].xout[-1]
                x = _x
            else:
                x = _new_x(res.xout, x, next_autonomous)
                y0 = res.yout[-1, :]
                if next_autonomous:
                    glob_x += res.xout[-1]
    if multimode:  # don't return defaultdict
        tot_nfo = [dict(nsys=oi+1, **_nfo) for _nfo in tot_nfo]
        return [Result(tot_x[idx], tot_y[idx], res[idx].params, tot_nfo[idx], odes[0])
                for idx in range(len(res))]
    else:
        tot_nfo = dict(nsys=oi+1, **tot_nfo)
        return Result(tot_x, tot_y, res.params, tot_nfo, odes[0])