It's certainly possible to turn these three second-order differential equations into six first-order ones. Your differentiation function will therefore have six inputs --- three coordinates and their three first derivatives --- and six outputs. (I think Matlab uses a vector of length six.)
The first three are trivial: the velocities your function spits out should be the velocities you put in:
\begin{align}
\frac d{dt} r &= r' \tag1\\
\frac d{dt} \theta &= \theta' \tag2\\
\frac d{dt} y &= y' \tag3
\end{align}
The other three derivatives you have to isolate from your equations of motion. The simplest route seems to want two auxiliary definitions. From your first equation,
$$
\Delta = r'' - \frac Mm y'' = r\theta'^2 + g\left(\sin\theta - \frac Mm \right).
$$
From (the derivative of) your second,
$$
\Sigma = r'' + y'' = -\frac d{dt}(r\theta') = - r'\theta' - r\theta''.
$$
(If I'm misreading your second equation and you mean a constant $R$ rather than the dynamical $r$, that's an easy substitution here, with $R' = 0$; note also $\theta''$ below.)
Then the derivatives of the velocities that you input are
\begin{align}
\frac d{dt} r' &= \frac{M\Sigma + m\Delta}{M+m} \tag4 \\
\frac d{dt} \theta' &= \frac{-2 r'\theta'- g\cos\theta}r \tag5 \\
\frac d{dt} y' &= \frac{M\Sigma - m\Delta}{2M} \tag6 \\
\end{align}
player100 points out in a comment that this system of equations actually isn't sensitive to a constant offset in $y$: the system is first-order in $y'$ rather than second-order in $y$. However, if $y$ is what interests you, then it's probably better to include it in your integrator than to leave it out. Integrating algorithms which make decisions about dynamically changing step sizes do so based on the error terms in the data that they have to look at.
You ask in a comment about when, in general, we can determine whether a system of equations is under/overdetermined; I think that's an interesting question, but too broad for this answer.
Note that a previous version of this answer put the "trivial" derivatives in the wrong place. Sorry about that --- it's been a long time since I ode45'd anything, and I'm rusty.
The equation can be rewritten as $(y')^2+yy'-x^2-x=0.$ One can solve $y'$ in terms of $y$ and $x.$ Notice that the above equation is equivalent to $(2y')^2+2y(2y')-4x^2-4x=0=(2y')^2+2y(2y')+y^2-(y^2+4x^2+4x)=(2y'+y)^2-(y^2+4x^2+4x),$ hence $(2y'+y)^2=y^2+4x^2+4x,$ implying $y^2\geq-(4x^2+4x).$ Notice that $4x^2+4x=4x^2+4x+1-1=(2x+1)^2-1,$ so $y^2\geq1-(2x+1)^2.$ This allows us to say that $2y'+y=\pm\sqrt{y^2+(2x+1)^2-1}.$
Best Answer
The goal here is to reduce the equation to something only in terms of $y,p,dp/dy$. If you can get $x$ to one side, it will go away after differentiation.
If you solve for $x$
$$ x = \frac{y - y^2p^3}{2p} = \frac{y}{2p}-\frac{y^2p^2}{2} $$
Then differentiate
$$ \frac{1}{p} = \frac{1}{2p} - \frac{y}{2p^2}\frac{dp}{dy} - yp^2 - y^2p\frac{dp}{dy} $$
Multiplying through by $2p^2$
$$ p + (y + 2y^2p^3)\frac{dp}{dy} + 2yp^4 = 0 $$
Then factor
$$ (1 + 2yp^3)\left(p + y\frac{dp}{dy}\right) = 0 $$
The equation is now solvable for $p = f(y)$