Lagrangian mechanics formulation of a simple free motion of two masses in uniform gravity field

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{dvec{r}_M}{dt}right|^2 , + , frac{m}{2},left|frac{dvec{r}_m}{dt}right|^2 , - , Mg,big(vec{r}_M cdot vec{e}_ybig) , - , mg,big(vec{r}_m cdot vec{e}_ybig)$$ with holonomic constriant $$big|vec{r}_M , - , vec{r}_mbig|^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{dvec{r}_G}{dt} , + , frac{m}{M+m} , frac{dvec{r}_l}{dt}right|^2 , + , frac{m}{2},left|frac{dvec{r}_G}{dt} , - , frac{M}{M+m} , frac{dvec{r}_l}{dt}right|^2 &, - ,g,Big(big(M,vec{r}_M , + , m , vec{r}_m big)cdot vec{e}_yBig) , =&, frac{M}{2} left|frac{dvec{r}_G}{dt}right|^2 , + , frac{M,m}{(M+m)}left(frac{dvec{r}_G}{dt}cdotfrac{dvec{r}_l}{dt}right) , + , frac{M,m^2}{2(M+m)^2}, left|frac{dvec{r}_l}{dt}right|^2 &, + , frac{m}{2} left|frac{dvec{r}_G}{dt}right|^2 , - , frac{M,m}{(M+m)}left(frac{dvec{r}_G}{dt}cdotfrac{dvec{r}_l}{dt}right) , + , frac{M^2,m}{2(M+m)^2} left|frac{dvec{r}_l}{dt}right|^2 &, - ,(M+m),g,Big(vec{r}_G cdot vec{e}_yBig) end{align} and after some regrouping and simplification the Lagrangian and the holonomic constraint become: $$L , = , frac{M,m}{2(M+m)}, left|frac{dvec{r}_l}{dt}right|^2 , + , frac{(M+m)}{2}, left|frac{dvec{r}_G}{dt}right|^2 , - , (M+m),g,Big(vec{r}_G cdot vec{e}_yBig) $$ $$big|vec{r}_lbig|^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}_mbig|^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^2vec{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}_ybig) = vec{0}times big(-(M+m), g , vec{e}_ybig) = vec{0}$$ i.e. the torque is zero, the equations of motion become begin{align} &frac{d^2vec{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^2vec{r}_G}{dt^2} , = , - , g, vec{e}_y & ,I , frac{domega}{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}_mbig|^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.