how to solve the system of differential equations for this particle? I'm trying to solve this problem
A particle of mass m moves under the action of gravity on the inner surface of a paraboloid of revolution $x^2+y^2=az$ which assumed frictionless. Obtain the equations of motion.
The Lagrangian in polar coordinates, assuming gravity towards the negative z axis, is $$L=\frac 1 2m\left(\dot\rho^2+\rho^2\dot\varphi^2+\dot z^2\right)-mgz;\qquad \ddot q\equiv\frac {d\dot q} {dt}\equiv\frac{d^2 q}{dt^2}$$constraints imposed condition implies $$f=\rho^2-az=0$$ Lagrange equations for this system are$$\frac d{dt}\left(\frac{\partial L}{\partial \dot q_i} \right)-\frac{\partial L}{\partial q_i}=\lambda\frac{\partial f}{\partial q_i}$$
$\lambda$ is a Lagrange multiplier. Then to $q_1=\rho$, $q_2=z$, $q_3=\varphi$ 
$$\fbox{$m\ddot\rho-m\rho\dot \varphi^2=2\rho\lambda$}$$
$$\fbox{$m\ddot z+mg=-a\lambda $}$$
$$\fbox{$ \frac{d}{dt}\left(\rho^2\dot\varphi\right)=0$}$$
$$\fbox{$\rho^2-az=0$}$$
$$\rho\in[0,+\infty),\quad\varphi\in[0,2\pi),\quad z\in(-\infty,+\infty)$$
I do not know how to solve this system of equations, and most important is to determine $\lambda$.
From system of equations I could deduce that
$$\rho^2\dot\varphi=c_0$$
$$m\ddot\rho-m\frac{c_0^2}{\rho^3}=2\rho\lambda$$
$$az\dot\varphi=c_0$$
$$\dot z\dot\varphi+z\ddot\varphi=0$$

in Cartesian coordinates
$$\fbox{$m\ddot y=2y\lambda$}$$
$$\fbox{$m\ddot x=2x\lambda$}$$
$$\fbox{$m\ddot z+mg=-a\lambda $}$$
$$\fbox{$x^2+y^2-az=0$}$$
 A: You seem to have successfully found a set of equations of motion for the particle. Yes you could eliminate some variables, but unless the problem says to do that, there's really no need.
If you want to solve this, well, things get a little more complex. We know this is an energy conserving and angular momentum conserving system, so we know that the equations of motion can be simplified to
$$
\rho^2\dot{\varphi} = c_0 \\
\frac{m}{2}\left(\dot{\rho}^2 + \rho^2\dot{\varphi}^2 + \dot{z}^2\right) + mgz = E\\
\rho^2 - az = 0.
$$
To derive this from the Lagrangian equations of motion, multiply the $\rho$ equation by $\dot{\rho}$, the $z$ equation by $\dot{z}$, and add them together. It will all simplify out to the conservation of energy given above.
Substituting in for $\dot{\varphi}$ and $z$ gives
$$
\frac{m}{2}\left[\left(1 + \frac{4\rho^2}{a^2}\right)\dot{\rho}^2 + \frac{c_0^2}{\rho^2}\right] + \frac{mg}{a}\rho^2 = E.
$$
This is a separable first order equation, so it's in principle solvable. We can get rid of a lot of the constants using $\rho(t) = (a/2)r(t\sqrt{2g/a})$:
$$
(1 + r^2)\dot{r}^2 + \frac{k}{r^2} + r^2 = \epsilon\;\;\;\;;\;\;\;\; k\equiv \frac{8c_0^2}{a^3 g}\,,\,\, \epsilon\equiv \frac{4E}{mga}
$$
While noticeably cleaner than before, it's still kind of a hot mess. Such is life in dynamics. Separating and using a $u$-substitution gives
$$
\frac{1}{2}\int_{r_0^2}^{r^2}\sqrt{\frac{1+u}{\epsilon u - u^2 - k}}du = t - t_0
$$
which is a solution, albeit an implicit one that may be hard to work with.
A: You usually don't find $\lambda$, but use your constraint to eliminate the unphysical degree of freedom (3 variables in your parametrization, 2 degrees of freedom in the system).
In more detail, differentiate the constraint twice, plug it in the $z$ equation of motion, multiply appropriately to equate the $\rho$ equation. You get 
$$
m  \ddot\rho -m \rho \dot\phi^2 = -2\rho\frac m a \left(g+ \frac 2 a(\dot\rho^2 + \rho\ddot\rho)\right)
$$
which leads to
$$
\ddot\rho-\rho\dot\phi^2 + 2\frac g a \rho+4\frac{\rho}{a^2}\left(\dot\rho^2 + \rho\ddot\rho\right)
=0.
$$
This is the equation of motion of the Lagrangian (eliminating $z$ from the start, $m=1$) 
$$
L=\frac 1 2 (\dot\rho^2 + \rho^2\dot\phi^2 +\frac{ 4}{ a^2} (\rho\dot\rho)^2)
-\frac g a \rho^2
$$
as it should be.
As to the solution of the equations, I do not know.
