How to find the equations of motion with a constraint? I am trying to find the equation of motion for a particle constrained to move on the surface defined by $S:z=\cos x+\sin y$ under the influence of gravity. I am working in the Cartesian Coordinate system, so the force of gravity is given by $F_g=-mg\hat{k}$ where $\hat{k}$ is a unit vector in the upwards z-direction.
Thus, the Lagrangian for this particle would be (I'm not sure if this is correct):
$$\mathcal{L}=\frac{1}{2}m\left(\dot{x}^2+\dot{y}^2+\dot{z}^2\right)-mgz+\lambda\left(\cos x+\sin y-z\right)$$
I am using the method of lagrange multipliers to solve for the equations of motion. Is there an easier way?
I then used the Euler-Lagrange equations to find my equations of motion:
$$\begin{align} m\ddot{x}&=-\lambda\sin x\\
m\ddot{y}&=\lambda \sin y\\
m\ddot{z}&=-mg-\lambda\\
z &=\cos x +\sin y
\end{align}$$
Now, I am not really sure what I need to do with the $\lambda$. I solved for it ($\lambda=-m(g+\ddot{z})$, and substituted it in the equations, leading me to these:
$$ \begin{align}
\ddot{x}&=(g+\ddot{z})\sin x \\
\ddot{y}&=(g+\ddot{z})\sin y \\
m\ddot{z}&=-mg-\lambda
\end{align} $$
But I still haven't completely isolated $\lambda$.
I then differentiated the constraint with respect to time twice, 
$$\ddot{z}=-\left(\dot{x}^2\cos x +\dot{y}^2\sin y\right)$$
And I substituted for $\ddot{z}$ in the equation solved for $\lambda$:
$$\lambda=-m\left(g-\left(\dot{x}^2\cos x +\dot{y}^2\sin y\right)\right)$$
I then eliminated $\lambda$ in the last equation:
$$m\ddot{z}=-mg+m\left(g-\left(\dot{x}^2\cos x +\dot{y}^2\sin y\right)\right)$$
$$\ddot{z}=-g+\left(g-\left(\dot{x}^2\cos x +\dot{y}^2\sin y\right)\right)$$
$$\ddot{z}=-\left(\dot{x}^2\cos x +\dot{y}^2\sin y\right)$$
Here are my questions: Is all of this correct?
Do you think there are analytical solutions to these differential equations?
 A: Deriving twice concerning $t$ the smooth constraint gives
$$
-x''(t) \sin (x(t))-x'(t)^2 \cos (x(t))+y''(t) \cos (y(t))-y'(t)^2 \sin (y(t))-z''(t)=0
$$
which jointly with
$$
\left\{
\begin{array}{rcl}
 -\lambda  \sin (x(t))-m x''(t)=0 \\
 \lambda  \cos (y(t))-m y''(t)=0 \\
 -\lambda -g m-m z''(t)=0 \\
\end{array}
\right.
$$
can be solved to $x''(t),y''(t),z''(t),\lambda$ giving
$$
\left\{
\begin{array}{rcl}
 x''(t)& = & \frac{\sin (x(t)) \left(-\cos (x(t)) x'(t)^2-\sin (y(t)) y'(t)^2+g\right)}{\cos ^2(y(t))+\sin ^2(x(t))+1} \\
 y''(t)& =& \frac{\cos (y(t)) \left(\cos (x(t)) x'(t)^2+\sin (y(t)) y'(t)^2-g\right)}{\cos ^2(y(t))+\sin ^2(x(t))+1} \\
 z''(t)& =& \frac{2 \left(\cos (x(t)) x'(t)^2+\sin (y(t)) y'(t)^2-g\right)}{\cos (2 x(t))-\cos (2 y(t))-4}-g \\
 \lambda & =& \frac{m \left(\cos (x(t)) x'(t)^2+\sin (y(t)) y'(t)^2-g\right)}{\cos ^2(y(t))+\sin ^2(x(t))+1} \\
\end{array}
\right.
$$
Attached the movement for initial conditions $x'(0) = y'(0)=0.01, x(0) = 0,y(0)=\pi/2$

A: To expand the commend of Valter Moretti, I would also use a set of coordinates which already satisfy the spatial constraints.
So I would start by parameterizing the position vector as
$$
\mathbf r(t) = \begin{pmatrix} x(t) \\ y(t) \\ z(t) \end{pmatrix}
= \begin{pmatrix} x(t) \\ y(t) \\ \lambda \cos(\mu\, x(t)) + \lambda\sin(\mu\, y(t)) \end{pmatrix},
$$
where I introduced the real parameters $\mu$, where $[\mu] = m^{-1}$, and $\lambda$, where $[\lambda] = m$, to take care of the units. By differentiation we get the velocity vector
$$
\dot{\mathbf{r}}(t) = \begin{pmatrix} \dot x(t) \\ \dot y(t) \\ -\lambda\mu \sin(\mu\, x(t))\ \dot x(t) + \lambda\mu\cos(\mu\, y(t))\ \dot y (t) \end{pmatrix}.
$$
The Lagrangian is then given by
$$
\begin{align}
L &= T(\dot{\mathbf{r}}) - U(\mathbf r)\\
&=\frac m 2 \dot{\mathbf{r}}^2 - mg\,\mathbf r\cdot \mathbf e _z \\
&=\frac m 2 (\dot x ^2 + \dot y ^2 + \lambda^2\mu^2(\sin^2(\mu x)\dot x^2 - 2\sin (\mu x)\cos(\mu y)\dot x \dot y +\cos^2(\mu x)\dot y^2)\\&\phantom{==} -mg\lambda(\cos(\mu x) + \sin(\mu y ))
\end{align}
$$
and the equations of motion are then given by
$$
\frac{\text{d}}{\text{d}t} \frac{\partial L}{\partial \dot x} - \frac{\partial L}{\partial x} = 0 \quad\text{and}\quad
\frac{\text{d}}{\text{d}t} \frac{\partial L}{\partial \dot y} - \frac{\partial L}{\partial y} = 0
$$
The function $z(t)$ is then found by using the solutions $x(t), y(t)$ in $z(t) = \lambda \cos(\mu\, x(t)) + \lambda\sin(\mu\, y(t))$.
