Unable to implement boundary conditions for a spherical wave equation coupled to a Laplace equation

Question

I am having trouble setting up in Mathematica v 12.1 the following system of coupled pdes on $r\in[0,10]$ $$ \left(\partial_t^2 - \frac{1}{r^2}\partial_r r^2 \partial_r \right) \theta(t,r) = -(1+2\Phi(t,r))\theta(t,r)$$ $$\frac{1}{r^2}\partial_r r^2 \partial_r \Phi(t,r) = \theta(t,r)^2 $$ For $\theta$ I have the initial conditions $$\theta(0,r)=e^{-r^2}, \qquad \partial_t\theta(0,r)=0.$$ and smoothness bc at origin and absorbing bc at $r=10$. $$\partial_r \theta(t,0)=0, \qquad \left.\partial_t(r \theta(t,r))\right|_{r=10}+\left.\partial_r(r \theta(t,r))\right|_{r=10}=0$$ The latter ensures no information comes in from outside the domain. $\Phi$ can be obtained by integrating the second equation so I don't believe I need to provide initial conditions for it. There are in principle two integration constants however so I must therefore impose two conditions, e.g. $$\partial_r \Phi(t,0)=0, \qquad \Phi(t,0)=-4$$

Below is how I've tried to implement it in Mathematica. (As far as I understand, the smoothness bc at the origin are taken care of automatically). I get the following error

tend = 100.;
rstart = 1/10^100;
rend = 10;
ic = {[Theta][0, r] == Exp[-(r)^2], (D[[Theta][t, r], t] /. t -> 0) == 0};
bc = {[CapitalPhi][t, 0] == -4};
{teta, phi} = NDSolveValue[Flatten@{-((2 r !(*SuperscriptBox[([Theta]),TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}],")"}],Derivative],MultilineFunction->None])[t, r] + r^2!(*SuperscriptBox[([Theta]), TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None])[t, r])/r^2) +!(*SuperscriptBox[([Theta]), TagBox[RowBox[{"(", RowBox[{"2", ",", "0"}], ")"}],Derivative],MultilineFunction->None])[t,r] == -(1 + 2 [CapitalPhi][t, r]) [Theta][t, r] + NeumannValue[-Derivative[1, 0][[Theta]][t, r] - [Theta][t, r]/r, r == rend] , ((2 r !(*SuperscriptBox[([CapitalPhi]), TagBox[RowBox[{"(",RowBox[{"0",",", "1"}], ")"}],Derivative],MultilineFunction->None])[t, r] + r^2 !(*SuperscriptBox[([CapitalPhi]), TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}],")"}],Derivative],MultilineFunction->None])[t, r])/r^2) == [Theta][t, r]^2, ic,bc}, {[Theta], [CapitalPhi]}, {t, 0, tend}, {r, rstart, rend}]
Animate[Plot[{teta[t, r], phi[t, r]}, {r, 0, rend},PlotRange -> {{0, rend}, {-1, 1}}], {t, 0, tend}]