'Find real solutions of the the equation ``pp(x) == y``.
Parameters
y : float, optional
Right-hand side. Default is zero.
discontinuity : bool, optional
Whether to report sign changes across discontinuities at
breakpoints as roots.
extrapolate : {bool, \'periodic\', None}, optional
If bool, determines whether to return roots from the polynomial
extrapolated based on first and last intervals, \'periodic\' works
the same as False. If None (default), use `self.extrapolate`.
Returns
roots : ndarray
Roots of the polynomial(s).
If the PPoly object describes multiple polynomials, the
return value is an object array whose each element is an
ndarray containing the roots.
Notes
This routine works only on real-valued polynomials.
If the piecewise polynomial contains sections that are
identically zero, the root list will contain the start point
of the corresponding interval, followed by a ``nan`` value.
If the polynomial is discontinuous across a breakpoint, and
there is a sign change across the breakpoint, this is reported
if the `discont` parameter is True.
Examples
Finding roots of ``[x**2 - 1, (x - 1)**2]`` defined on intervals
``[-2, 1], [1, 2]``:
>>> from scipy.interpolate import PPoly
>>> pp = PPoly(np.array([[1, -4, 3], [1, 0, 0]]).T, [-2, 1, 2])
>>> pp.roots()
array([-1.,  1.])'

def solve(self, y=0.0, discontinuity=True, extrapolate=None):

    if (extrapolate is None):
        extrapolate = self.extrapolate
    self._ensure_c_contiguous()
    if np.issubdtype(self.c.dtype, np.complexfloating):
        raise ValueError('Root    finding    is    only    for    real-valued    polynomials')
    y = float(y)
    r = _ppoly.real_roots(self.c.reshape(self.c.shape[0], self.c.shape[1], (-1)), self.x, y, bool(discontinuity), bool(extrapolate))
    if (self.c.ndim == 2):
        return r[0]
    else:
        r2 = np.empty(prod(self.c.shape[2:]), dtype=object)
        for (ii, root) in enumerate(r):
            r2[ii] = root
        return r2.reshape(self.c.shape[2:])
