Is this differential equation solvable in MATLAB? \begin{align}
    r'' - r θ'^2 &= g \sinθ - (M/m)(g - y'')
\\    
    r' + Rθ' + y' &= 0
  \\  
    rθ'' + 2r'θ' + g \cosθ &= 0
\end{align}
Sorry may I know how I can solve this with ode45? I am attempting to convert all three equations into 1st Order Coupled equations, but I can't seem to get the second equation into the correct substitution. Is it possible? Or do I have to use another ODE solver (which?)
Edit: I attempted rob's solution, and have gotten the following code
Taking M=3m, g=9.81, R=4
function dydt = odefun(t,y)
dydt = zeros(9,1)
dydt(1) = y(7)
dydt(2) = y(8)
dydt(3) = y(9)
dydt(4) = sqrt((y(8)-9.81*sin(y(4))-3*9.81+3*y(9))/(2*y(1)))
dydt(5) = -(9.81*cos(y(4))+y(5)*y(7))/(2*y(1))
dydt(6) = -y(2)-y(1)*4
end

but received the following error message
Not enough input arguments.
Error in <filename> (line 3)
dydt(1) = y(7)

It seems matlab does not accept the code if I do not define (θ'')' = ___. Sorry if this is a dumb misunderstanding, I am still trying to figure out how MATLAB works (any resources to recommend?)
I have also attempted 
function dydt = odefun(t,y)
dydt = zeros(9,1)
dydt(1) = y(7)
dydt(2) = y(8)
dydt(3) = y(9)
dydt(4) = sqrt((y(8)-9.81*sin(y(4))-3*9.81+3*y(9))/(2*y(1)))
dydt(5) = -(9.81*cos(y(4))+y(5)*y(7))/(2*y(1))
dydt(6) = -y(2)-y(1)*4 % radius here!
dydt(4) = y(1) % are these 3 more lines correct logically, they seem to 
dydt(5) = y(2) % let matlab run, but for now I am just getting constants
dydt(6) = y(3) % for all solutions
end

is this the right code?
 A: 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.
A: Take $r-y = u$ and $r+y = v$ then
\begin{align}
 u''+\frac{(v+u)\theta'^2}{2} &= g \sin\theta -g  
\\ v'+ \theta' &= 0 
\\ \frac{(v+u)}{2} \theta''+ \frac{(v'+u')}{2}+g \cos\theta &= 0
\end{align}
From the second equation, $ \theta = c-v$ for some constant $c$. Then take $x_1 =u, x_2=v, x_3= u', x_4=v' $ so for numerical solution
\begin{align}
   \frac{dx_1}{dt} &= x_3 \\
   \frac{dx_2}{dt} &= x_4 \\
   \frac{dx_3}{dt} &= -\frac{(x_1+x_2)x_4^2}{2}+g\sin(c-x_2)-g \\
   \frac{dx_4}{dt} &= \frac{x_3+x_4}{x_1+x_2}+\frac{2g\cos(c-x_2)}{x_1+x_2} 
\end{align}
For initial conditions you need only five of them: four for differential equations and one for determining $c$. Also there seems to be a potential annoyance around $x_1=-x_2$, i.e. $r=0$. 
