Runge-Kutta on coupled equations

Question

I'm still pretty new to Mathematica and I'm having a few issues solving these coupled equations. I've attempted using NDSolve with method "ExplicitRungeKutta" as well as with no method, and have had zero success so far. I came across a few different answers on here (e.g. Solving a system of ODEs with the Runge-Kutta method) that I attempted to implement, but haven't had any luck yet. Any suggestions would be greatly appreciated.

Thanks

ClearAll["Global'*"]
ClearAll["@"]
ClearAll[\[Phi], Tru, Trp, ruf, ruv, ru, rp, \[CapitalOmega], \
\[CapitalOmega]0, \[Rho]su, rb0, H, h, \[Sigma], cpb, \[Rho]0, Ru, \
Rmix, \[Eta], wh2, r, r\[Prime], P, \[Upsilon], Qfis, step, keff, \
whl, whv]
Constants -> {Tru, Trp, 
  ru, \[CapitalOmega], \[CapitalOmega]0, \[Rho]su, H, h, \[Sigma], 
  cpb, \[Eta], wh2, P}

step = 0.001;
Qfis = 200;(* MeV/fission*)
\[Upsilon] = 2.464; (*average number of neutrons produced per fission*)


P = 10.0*10^6; (*Watts*)
keff = 1.00598; Mu0 = 5.99685*^4; Tru = 800; Tboil = 4404; Tmelt = \
1405.3; Trp = 5500; ru = 0.05;
rp = 0.03; PB0 = 30*^6; \[Sigma] = 1500*^-3; rb0 = 100*^-6; Rsu = 
 8.314/238;
Rsh2 = 8.314/2.016; ku = 27; H = 1; \[Eta] = 0.5; \[Rho]h2 = 70.85;
rlist = Range[rp, ru, step];
KCODEnorm = (\[Upsilon]*P)/((1.602*10^-13)*Qfis*keff);
Edep[r_] = KCODEnorm*(1.15455*^-3);
\[Rho]0[r_] = Mu0/(\[Pi] (ru^2 - rp^2) H);(*g m^-3*)
\[Phi][r_] = (Edep[r]*\[Rho]0[r]*(1.602*10^-13))/
  ku; (*1 MeV = 1.602E-13 W*s *)
ansol = FullSimplify[
   DSolve[{\[Phi][r] + Derivative[1][y][r]/
       r + (y^\[Prime]\[Prime])[r] == 0, y[rp] == Trp, y[ru] == Tru}, 
    y[r], r]];
Tu[r_] = FullSimplify[Evaluate[y[r] /. ansol], 
  ru > 0 && rp > 0 && r > 0]

rsolV = Solve[Tu[r] == Tboil, r]
rsolF = Solve[Tu[r] == Tmelt, r]
ruv = r /. rsolV[[2]];
ruf = r /. rsolF[[2]];

\[Rho]su = 19100;
\[Rho]l[r_] = (17270 - 1.4485 (Tu[r] - 1408))
Pl[r_, \[CapitalOmega]_] = 
  Integrate[\[Rho]l[
     r\[Prime]] r\[Prime] \[CapitalOmega]^2, {r\[Prime], 0, r}];
\[Rho]v[r_, \[CapitalOmega]_] = (Pl[ruv, \[CapitalOmega]]/(
   Rsu*Tu[ruv])) Exp[(\[CapitalOmega]^2/(2*Rsu)) (r^2/Tu[r] - ruv^2/
      Tu[ruv] )]

Pv[r_, \[CapitalOmega]_] = \[Rho]v[r, \[CapitalOmega]][[1]]*Rsu*Tu[r];
\[Rho][r_, \[CapitalOmega]_] = 
  Piecewise[{{\[Rho]v[r, \[CapitalOmega]][[1]], 
     rp <= r < ruv}, {\[Rho]l[r][[1]], ruv <= r < ruf}, {\[Rho]su, 
     ruf <= r <= ru}}];
Pu[r_, \[CapitalOmega]_] = 
  Piecewise[{{Pv[r, \[CapitalOmega]][[1]], 
     rp <= r <= ruv}, {Pl[r, \[CapitalOmega]][[1]], 
     ruv <= r <= ruf}, {PB0, ruf <= r <= ru}}];

\[CapitalOmega]sol = 
 Solve[Pl[ruf, \[CapitalOmega]] == 
   PB0 - (2 \[Sigma])/rb0, \[CapitalOmega]]
\[CapitalOmega]0 = \[CapitalOmega] /. \[CapitalOmega]sol[[2]];

(**********************************************
Solve coupled equations
********************************************)

cpb = 15(*J/gK*); h = 1;
Mb = PB0/(Rsh2 Tu[ruf]) (4/3 \[Pi]) rb0^3;
C0 = (3 Mb Rsh2)/(4 \[Pi]);

func = {
  Tb'[t] == (((4*\[Pi]*(rb[t])^2*h)/
        cpb )*(Tu[Rb[t]] - Tb[t]) - (\[Rho]b'[t])*Tb[t])/\[Rho]b[t],
  Rb''[t] == -Mb^-1*
    Integrate[(4 \[Pi] r\[Prime]^2 (Dt[Pu[r, \[CapitalOmega]0], r] + 
           Pu[r, \[CapitalOmega]0]/r) /. (r -> Rb[t] + r\[Prime])) // 
      FullSimplify, {r\[Prime], 0, rb[t]}, 
     Assumptions -> Element[{Rb[t], rb[t]}, Reals]],

  rb[t] == (2 \[Sigma])/((C0*Tb[t])/(rb[t])^3 - 
     Pu[Rb[t], \[CapitalOmega]0]),
  \[Rho]b[t] == C0/(Rsh2*(rb[t])^3)
  }

yinit = {Tb[0] == Tu[ruf], Rb[0] == ruf, rb[0] == rb0, Rb'[0] == 0, 
   rb'[0] == 0, Rb''[0] == 0};
y = {Tb[t], Rb[t], rb[t], \[Rho]b[t]};

tf = 1*^-3;
tstep = 1*^-6;
(*sol=NDSolve[Join[func,yinit],y,{t,0,tf},Method\[Rule]\
"ExplicitRungeKutta", "StartingStepSize"\[Rule]tstep]*)
NDSolve[Join[func, yinit], y, {t, 0, tf}]