Lagrangian mechanics formulation of a simple free motion of two masses in uniform gravity field As a part of larger project, I decided to test my Lagrangian formulation of simple system of two rigidly connected point masses as indicated below.

I introduce the generalized coordinates vector $\bar{q} = [x, y, q]^T$, where $x,y$ are coordinates of position of mass $M$ and $q$ is the angular deflection from horizontal of the massless rod of length $l$, in the middle of which mass $m$ is attached. Then, the position coordinates of $m$ are given by
$$ x_m = x + \dfrac{l}{2}\cos(q) $$
$$ y_m = y + \dfrac{l}{2}\sin(q) $$
Having this, I defined kinetic ($T$) and potential ($V$) energies
$$ T = \dfrac{1}{2}M(\dot{x}^2 + \dot{y}^2) +  \dfrac{1}{2}m(\dot{x}_m^2 + \dot{y}_m^2)$$
$$ V = g \left( My + m\left(y+\dfrac{l}{2}\sin(q)\right)  \right) $$
Then, the Lagrangian $L$ and mechanical energy of the system $E$ are
$$ L = T-V$$
$$ E = T+V$$
Further, the derivations were done by my by hand (a few times to double check) and using SymPy (symbolic math package of Python) and they match. Upon integrating the EOM, I receive a nicely looking plots for state vector variables $\bar{Q} = [x,y,q,\dot{x},\dot{y},\dot{q}]^T$, if my initial conditions involve null $q$ nd $\dot{q}$. However, if they are non-zero, the motion looks highly non-physical visually (I created a Matlab animation of the motion) and additionally the mechanical energy of the system is not adding up to be constant (see pic below for $\bar{Q}_0 = [1,1000,1,1,0,-0.1]^T$).

My question is therefore, are my assumptions and initial formulation correct for the given situation? I am quite new to Lagrangian formalism and afters hours of tackling this seemingly simple problem, started to think there might be some rule of Lagrangian mechanics I violate. All classical mechanics problems I was able to find online are not free motion of multibody systems, therefore I was not able to check with actual solution for problem like this.
To give some more context, my intention is to expand this problem to put the multibody system in orbital flight (central point-mass gravity field) and increase complexity of structure, leading to n-link pendulum attached to bus mass $M$. If there is a fundamental error in the proceedings I presented above here, I assume I would not be able to extend the formulation to more complex system I just described.
I will be eternally grateful for any help / advice! I will be also happy to provide more details on my solution.
EDIT
To clarify, I post my EL equations:
$$ M \frac{d^{2}}{d t^{2}} x{\left(t \right)} - 0.5 m \left(l \sin{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} q{\left(t \right)} + l \cos{\left(q{\left(t \right)} \right)} \left(\frac{d}{d t} q{\left(t \right)}\right)^{2} - 2 \frac{d^{2}}{d t^{2}} x{\left(t \right)}\right)= 0 $$
$$ M \frac{d^{2}}{d t^{2}} y{\left(t \right)} + g \left(M + m\right) + 0.5 m \left(- l \sin{\left(q{\left(t \right)} \right)} \left(\frac{d}{d t} q{\left(t \right)}\right)^{2} + l \cos{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} q{\left(t \right)} + 2 \frac{d^{2}}{d t^{2}} y{\left(t \right)}\right) = 0 $$
$$ l m \left(0.5 g \cos{\left(q{\left(t \right)} \right)} + 0.25 l \frac{d^{2}}{d t^{2}} q{\left(t \right)} - 0.5 \sin{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} x{\left(t \right)} + 0.5 \cos{\left(q{\left(t \right)} \right)} \frac{d^{2}}{d t^{2}} y{\left(t \right)}\right) = 0$$
$$ \begin{bmatrix} M+m & 0 & -\dfrac{l}{2}m\sin(q) \\ 0 & M+m  & \dfrac{l}{2}m\cos(q) \\ -\dfrac{l}{2}m\sin(q) & \dfrac{l}{2}m\cos(q) & m\dfrac{l^2}{4}    \end{bmatrix} \begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \ddot{q}    \end{bmatrix} = \begin{bmatrix}  \dfrac{l}{2}m\cos(q)\dot{q}^2 \\ -g(M+m) + \dfrac{l}{2}m\sin(q)\dot{q}^2 \\ -\dfrac{l}{2}gm\cos(q)   \end{bmatrix} $$
 A: Remark. This case generalizes to your real goal to have a system of more than two mass-points attached to each other by rods where the rods can rotate relative to each other. Basically, in a uniform force field like the gravitational force field near the surface of the earth the problem becomes equivalent  to the geodesic motion of a point onto a torus of dimension equal to the number of rods.
Lagrangian approach. First, you place an inertial coordinate frame (say a frame attached to the ground) in the plane of the system. Denote it by $O\,\vec{e}_x\,\vec{e}_y$, where the coordinate vectors are of unit length and perpendicular to each other. Moreover, $O\,\vec{e}_x\,\vec{e}_y$ is oriented so that $\vec{e}_y$ is the vertical vector along which gravitational acceleration is $\vec{g} = - \, g \, \vec{e}_y$. The position vectors of the mass-points, originating from point $O$ and ending at the either of the mass points, are
$\vec{r}_M$ and $\vec{r}_m$. The Lagrangian is then
$$L\, =\, \frac{M}{2}\,\left|\frac{d\vec{r}_M}{dt}\right|^2 
\, + \,  \frac{m}{2}\,\left|\frac{d\vec{r}_m}{dt}\right|^2 
\, - \, Mg\,\big(\vec{r}_M \cdot \vec{e}_y\big)
\, - \, mg\,\big(\vec{r}_m \cdot \vec{e}_y\big)$$
with holonomic constriant
$$\big|\vec{r}_M \, - \, \vec{r}_m\big|^2 \, = \, \frac{l^2}{4}$$
Perform the linear chainge of coordinates:
\begin{align}
&\vec{r}_G \, =\, \frac{M}{M+m} \, \vec{r}_M \, + \,  \frac{m}{M+m} \, \vec{r}_m\\
&\vec{r}_l \, = \, \vec{r}_M \, - \, \vec{r}_m 
\end{align}
with inverse
\begin{align}
&\vec{r}_M \, =\, \vec{r}_G \, + \, \frac{m}{M+m} \, \vec{r}_l\\
&\vec{r}_m \, = \, \vec{r}_G \, - \, \frac{M}{M+m} \, \vec{r}_l 
\end{align}
In the new coordinates, the Lagrangian becomes
\begin{align}
L \, =& \, \frac{M}{2}\,\left|\frac{d\vec{r}_G}{dt} \, + \, \frac{m}{M+m} \, \frac{d\vec{r}_l}{dt}\right|^2 
\, + \,  \frac{m}{2}\,\left|\frac{d\vec{r}_G}{dt} \, - \, \frac{M}{M+m} \, \frac{d\vec{r}_l}{dt}\right|^2 \\
&\, - \,g\,\Big(\big(M\,\vec{r}_M \, + \, m \, \vec{r}_m \big)\cdot \vec{e}_y\Big) \\
\, =&\, \frac{M}{2} \left|\frac{d\vec{r}_G}{dt}\right|^2 \, + \, \frac{M\,m}{(M+m)}\left(\frac{d\vec{r}_G}{dt}\cdot\frac{d\vec{r}_l}{dt}\right) \, + \, \frac{M\,m^2}{2(M+m)^2}\, \left|\frac{d\vec{r}_l}{dt}\right|^2\\
&\, + \, \frac{m}{2} \left|\frac{d\vec{r}_G}{dt}\right|^2 \, - \, \frac{M\,m}{(M+m)}\left(\frac{d\vec{r}_G}{dt}\cdot\frac{d\vec{r}_l}{dt}\right) \, + \, \frac{M^2\,m}{2(M+m)^2} \left|\frac{d\vec{r}_l}{dt}\right|^2\\
&\, - \,(M+m)\,g\,\Big(\vec{r}_G \cdot \vec{e}_y\Big) 
\end{align}
and after some regrouping and simplification the Lagrangian and the holonomic constraint become:
$$L \, = \, \frac{M\,m}{2(M+m)}\, \left|\frac{d\vec{r}_l}{dt}\right|^2
\, + \, \frac{(M+m)}{2}\, \left|\frac{d\vec{r}_G}{dt}\right|^2 \, - \, (M+m)\,g\,\Big(\vec{r}_G \cdot \vec{e}_y\Big) $$
$$\big|\vec{r}_l\big|^2 = \frac{l^2}{4}$$
Introduce the generalized coordinates:
$$\vec{r}_l \, = \, \frac{l}{2}\, \cos(q) \, \vec{e}_x \, + \, \frac{l}{2}\, \sin(q) \, \vec{e}_y$$
$$\vec{r}_G \, = \, x_G \, \vec{e}_x \, + \, y_G\, \vec{e}_y$$ which parametrixe the holonimc constrints (they tirivialize it). Then, the Lagrnagian simplifies to
$$L \, = \, \frac{M\,m\ l^2}{8(M+m)}\, \left(\frac{dq}{dt}\right)^2
\, + \, \frac{(M+m)}{2}\, \left(\frac{dx_G}{dt}\right)^2 \, + \,  \frac{(M+m)}{2}\, \left(\frac{dy_G}{dt}\right)^2 \, - \, (M+m)\,g\,y_G $$
Consequently, the Euler-Lagrange equations are
\begin{align}
&\frac{d^2{x}_G}{dt^2} \, = \, 0 \\
&\frac{d^2{y}_G}{dt^2} \, = \, - \, g\,\\
& \frac{d^2q}{dt^2} \, = \, 0
\end{align}
You can solve them right away
\begin{align}
&x_G \, = \, x_G(0) \, + \, v_x\, t \\
&y_G \, = \, y_G(0) \, + \, v_y\, t \, - \, \frac{g}{2} \, t^2\,\\
&q \, = \, q(0) \, + \, \omega_0 \, t 
\end{align}
where $x_G(0), \, y_G(0), \, q(0)$ are constants describing the initial configuration of the rod and $v_x, \, v_y, \, \omega_0$ are constants describing the initial velocities of the rod.
The rigid body interpretation. The special case you are modelling in your OP can be thought of as a solid body, basically a rigid rod, which moves in a plane and whose mass is concentrated at its ends. First, you place an inertial coordinate frame (say a frame attached to the ground). Denote it by $O\,\vec{e}_x\,\vec{e}_y\,\vec{e}_z$, where the coordinate vectors are of unit length and perpendicular to each other. Moreover, $O\,\vec{e}_x\,\vec{e}_y\,\vec{e}_z$ is oriented so that $\vec{e}_y$ is the vertical vector along which gravitational acceleration is $\vec{g} = - \, g \, \vec{e}_y$, the system moves in the coordinate plane $O\,\vec{e}_x\,\vec{e}_y$, just like on your picture, and $\vec{e}_z$ is pointing perpendicularly to the picture.
Denote by $\vec{r}_M$ the position vector, pointing from the origin $O$ of the coordinate frame to the mass-point $M$ and by $\vec{r}_m$ the position vector, pointing from the origin $O$ of the coordinate frame to the mass-point $m$. The condition that the two mass points are connected by a mass-less rod of length $\frac{l}{2}$, can be written as the quadratic equation.
$$\big|\vec{r}_M - \vec{r}_m\big|^2 = \frac{l^2}{4}$$ Let $G$ be the center of gravity of the rod, i.e. the center of gravity, which is
$$\vec{r}_G = \frac{M}{M+m} \, \vec{r}_M \, + \,  \frac{m}{M+m} \, \vec{r}_m$$
The force that acts on the system of the two mass points on the rod is gravity and is applied to the center of mass $\vec{r}_G$ of the rod. The general equations of motion of a rigid body can be written as
\begin{align}
&\frac{d^2\vec{r}_G}{dt^2} \, = \, - \, g\, \vec{e}_y\\
& \frac{d}{dt}\,I \,\omega\, \vec{e}_z \, = \, \vec{T}
\end{align} where $\omega \, \vec{e}_z$ is the angular velocity of the rod and $$I = \frac{M\,m\,l^2}{4\,(M+m)}$$ is its moment of inertia along the $O\, \vec{e}_z$ axis with respect to the center of mass $G$ of the rod. The other moments of inertia are zero because this is a one dimensional rod. The vector $\vec{T}$ is the sum of the torques of all forces acting on the rod, calculated with respect to the center of mass $G$. However, there is only one force, the gravity force $-(M+m)\, g \, \vec{e}_y$, acting on the rod and it is applied to the center of gravity $G$. Since $$\vec{T} = \vec{GG} \times \big(-(M+m)\, g \, \vec{e}_y\big) = \vec{0}\times \big(-(M+m)\, g \, \vec{e}_y\big) = \vec{0}$$ i.e. the torque is zero, the equations of motion become
\begin{align}
&\frac{d^2\vec{r}_G}{dt^2} \, = \, - \, g\, \vec{e}_y\\
& \frac{d}{dt}\,I \,\omega\, \vec{e}_z \, = \, \vec{0}
\end{align}
As you can see, these equations are decoupled and can be written as
\begin{align}
&\frac{d^2\vec{r}_G}{dt^2} \, = \, - \, g\, \vec{e}_y\\
& \,I \, \frac{d\omega}{dt} \, = \, {0}
\end{align}
But, the angular velocity, in your notations, is simply $\omega = \frac{dq}{dt}$ and also
$$\vec{r_G} = x_G \, \vec{e}_x + y_G\, \vec{e}_y$$ so
\begin{align}
&\frac{d^2{x}_G}{dt^2} \, = \, 0 \\
&\frac{d^2{y}_G}{dt^2} \, = \, - \, g\,\\
& \,I \, \frac{d^2q}{dt^2} \, = \, {0}
\end{align}
You can solve them right away
\begin{align}
&x_G \, = \, x_G(0) \, + \, v_x\, t \\
&y_G \, = \, y_G(0) \, + \, v_y\, t \, - \, \frac{g}{2} \, t^2\,\\
&q \, = \, q(0) \, + \, \omega_0 \, t 
\end{align}
where $x_G(0), \, y_G(0), \, q(0)$ are constants describing the initial configuration of the rod and $v_x, \, v_y, \, \omega_0$ are constants describing the initial velocities of the rod.
If I have time, I can get to these equations from Lagrangian point of view too. Basically you just have to change coordinates of your Lagrangian and the holonomic constrant $\big|\vec{r}_M - \vec{r}_m\big|^2 = \frac{l^2}{4}$ like this:
\begin{align}
&\vec{r}_G \, =\, \frac{M}{M+m} \, \vec{r}_M \, + \,  \frac{m}{M+m} \, \vec{r}_m\\
&\vec{r}_l \, = \, \vec{r}_M \, - \, \vec{r}_m 
\end{align}
Remark. This case does not generalize to your real goal to have a system of more than two mass-points attached to each other by rods where the rods can rotate relative to each other.
A: This is my solution for your problem:
The Kinetic Energy is:
$$T=\frac{1}{2}\,M(\dot x^2+\dot y^2)+\frac 12\,m\,\left[\left(\dot x -\frac 12 l \sin(\varphi)\,\dot \varphi\right)^2
+\left(\dot y +\frac 12 \,l\cos(\varphi)\dot \varphi\right)^2\right]
$$
and
the potential energy
$$U=M\,g\,y+m\,g\,\left(y +\frac 12 l\sin(\varphi)\right)$$
the generalized coordinates are:
$$\vec q=\begin{bmatrix}
   x \\
   y \\
   \varphi \\
 \end{bmatrix}$$
The EOM's:
$$\ddot x-\frac 12 \frac{\cos(\varphi)\,m\,l}{M+m}\,\dot \varphi^2=0$$
$$\ddot y-\frac 12 \frac{\sin(\varphi)\,m\,l}{M+m}\,\dot \varphi^2+g=0$$
$$\ddot \varphi=0$$
the numerical simulation  give you for $E=T+U$

the initial conditions are all zero except for $~\varphi(0)=0.3$
thus: if you have this equations, the total energy $E=T+U~$ is constant as it should be.
edit
in case that $\varphi(0)=\varphi_0$ and all other initial conditions are zero, the solution of the equations of motion are:
$$x(\tau)=0~,y(\tau)=-\frac 12\,g\,\tau^2~,\varphi(\tau)=\varphi_0 $$
and
$$E=\frac 12 m\,g\,l\,\sin(\varphi_0)=~\text{constant}$$
The equations of motion:
from Euler Lagrange  you obtain this three equations
$$\left( M+m \right) { \ddot x}-1/2\,ml\sin \left( \varphi  \right) 
\ddot \varphi -1/2\,ml\cos \left( \varphi  \right) {\dot \varphi }^{2}
=0$$
$$\left( M+m \right) {\ddot y}+1/2\,ml\cos \left( \varphi  \right) 
\ddot\varphi -1/2\,ml\sin \left( \varphi  \right) {\dot \varphi }^{2}+M
g+mg=0
$$
and
$$-1/4\,ml \left( 2\,\sin \left( \varphi  \right) { \ddot x}-2\,\cos
 \left( \varphi  \right) {\ddot y}-\ddot\varphi \,l-2\,g\cos
 \left( \varphi  \right)  \right) 
=0$$
solve those equations for $\ddot x~,\ddot y~,\ddot \varphi$ you obtain my equations above
