This is the problem I solved sometime ago. Same problem as you show. The angle $\theta$ used is measured from x-axis, positive anti-clockwise, as it was simpler to do so, but it does not affect the solution ofcourse. Since there is no Mathematica stuff in this, I can add a Manipulate later on if needed to show the ball falling of?
There are two coordinates $r,\theta$ (polar) which is the position vector of the ball, and one constraint
\begin{equation}
f\left( r,\theta\right) =r-R=0 \tag{1}
\end{equation}
$R$ above is the radius of the hemisphere.
Now we set up the equations of motion for $m$
\begin{align*}
T & =\frac{1}{2}m\left( \dot{r}^{2}+r^{2}\dot{\theta}^{2}\right) \\
U & =mgr\sin\theta\\
L & =T-U\\
& =\frac{1}{2}m\left( \dot{r}^{2}+r^{2}\dot{\theta}^{2}\right)
-mgr\sin\theta
\end{align*}
Hence the Euler-Lagrangian equations are (we need to add contraint)
\begin{align}
\frac{d}{dt}\frac{\partial L}{\partial\dot{r}}-\frac{\partial L}{\partial
r}+\lambda\frac{\partial f}{\partial r} & =0\tag{2}\\
\frac{d}{dt}\frac{\partial L}{\partial\dot{\theta}}-\frac{\partial L}
{\partial\theta}+\lambda\frac{\partial f}{\partial\theta} & =0 \tag{3}
\end{align}
But
\begin{align*}
\frac{d}{dt}\frac{\partial L}{\partial\dot{r}} & =m\ddot{r}\\
\frac{\partial L}{\partial\dot{\theta}} & =mr^{2}\dot{\theta}\\
\frac{d}{dt}\left( \frac{\partial L}{\partial\dot{\theta}}\right) &
=m\left( 2r\dot{r}\dot{\theta}+r^{2}\ddot{\theta}\right) \\
\frac{\partial L}{\partial r} & =mr\dot{\theta}^{2}-mg\sin\theta\\
\frac{\partial L}{\partial\theta} & =-mgr\cos\theta\\
\frac{\partial f}{\partial r} & =1\\
\frac{\partial f}{\partial\theta} & =0
\end{align*}
Hence (2) becomes
\begin{equation}
m\ddot{r}-mr\dot{\theta}^{2}+mg\sin\theta+\lambda=0 \tag{4}
\end{equation}
And (3) becomes
\begin{align}
m\left( 2r\dot{r}\dot{\theta}+r^{2}\ddot{\theta}\right) +mgr\cos\theta &
=0\nonumber\\
r\ddot{\theta}+2\dot{r}\dot{\theta}+g\cos\theta & =0 \tag{5}
\end{align}
We now need to solve (1,4,5) for $\lambda$. Now we have to apply the constrain
that $r=R$ in the above to be able to solve (4,5) equations. Therefore, (4,5)
becomes
\begin{align}
-mR\dot{\theta}^{2}+mg\cos\theta+\lambda & =0\tag{4A}\\
R\ddot{\theta}+g\cos\theta & =0 \tag{5A}
\end{align}
Where (4A,5A) were obtained from (4,5) by replacing $r=R$ and
$\dot{r}=0$ and $\ddot{r}=0$ since we are using that $r=R$ which is constant
(the radius).
From (5A) we see that this can be integrated giving
\begin{equation}
R\dot{\theta}^{2}+2g\sin\theta+c=0 \tag{6}
\end{equation}
Where $c$ is constant. Since if we differentiate the above with time, we
obtain
\begin{align*}
2R\dot{\theta}\ddot{\theta}+2g\dot{\theta}\cos\theta & =0\\
R\ddot{\theta}+g\cos\theta & =0
\end{align*}
Which is the same as (5A). Therefore from (6) we find $\dot{\theta}^{2}$ to
use in (4A). Hence from (6)
$$
\dot{\theta}^{2}=-2\frac{g}{R}\sin\theta+c
$$
To find $c$ we use initial conditions. At $t=0$, $\theta=90^{0}$ and
$\dot{\theta}\left( 0\right) =0$ hence
$$
c=2\frac{g}{R}
$$
Therefore
\begin{align*}
\dot{\theta}^{2} & =-2\frac{g}{R}\sin\theta+2\frac{g}{R}\\
& =2\frac{g}{R}\left( 1-\sin\theta\right)
\end{align*}
Plugging the above into (4A) in order to find $\lambda$ gives
\begin{align*}
-mR\left( 2\frac{g}{R}\left( 1-\sin\theta\right) \right) +mg\sin
\theta+\lambda & =0\\
\lambda & =m\left( 2g\left( 1-\sin\theta\right) \right) -mg\sin\theta\\
\lambda & =2mg-2mg\sin\theta-mg\sin\theta\\
& =mg\left( 2-3\sin\theta\right)
\end{align*}
Now that we found $\lambda\,,$we can find the constraint force in the radial
direction
\begin{align*}
N & =\lambda\frac{\partial f}{\partial r}\\
& =mg\left( 2-3\sin\theta\right)
\end{align*}
The particle will leave when $N=0$ which will happen when
\begin{align*}
2-3\sin\theta & =0\\
\theta & =\sin^{-1}\left( \frac{2}{3}\right) \\
& =41.8^{0}
\end{align*}
Therefore, the angle from the vertical is
$$
90-41.8=48.2^{0}
$$
r*θ''[t] == g*Sin[θ[t]], the solution will be $\theta(t)=0$, so the output of the first code sample is correct. – xzczd Jun 06 '20 at 11:35L = T-V.. – A little mouse on the pampas Jun 06 '20 at 11:37θ[t0]==0, why doesn'tL = 1/2 m (R θ'[t0])^2+m*g*R (1 - Cos[θ[t]]) - m*g*R*Cos[θ[t]]work?", as someone never learned about Lagrangian mechanics, I don't know the answer, either. Once again, I believe this should be asked in physics.SE. – xzczd Jun 08 '20 at 06:17L=m*g*R (1 - Cos[θ[t]]) - m*g*R*Cos[θ[t]]for analysisI . I always feel insecure because of the lack of authoritative arguments. So I'm interested in the textbook's authoritative analysis of this method. There are too few people in the SE physics section to discuss problems actively. Thank you very much. That's it. – A little mouse on the pampas Jun 08 '20 at 22:24