Consider the following system
Here, two immovable objects with mass $M$ are positioned $2d$ distance apart (in an empty universe). Meanwhile, an object of mass $m$ is placed somewhere above them on the same plane.
I found this system interesting because at first it resembled a simple harmonic oscillator, but on further inspection since the force of gravity doesn't increase linearly, it is not quite one (not dissimilar to a pendulum).
I'd like to model the motion of $m$ using the angle $\theta$ it makes with one of the unmovable masses $M$, but using Newtonian Mechanics, this seems quite nuanced.
I know that if we let $r$ be the distance between $m$ and one of the unmovable masses $M$, the force on the object that drives it forwards (along the blue dotted line) is (the factor of 2 comes from the fact that there's two unmovable masses pulling on it) $$F=\frac{2mMG}{r^2}\cdot\sin(\theta)$$ Since the blue distance squared plus $d^2$ is $r^2$, and $d\tan(\theta)$ gives the blue distance, we an rewrite the previous force expression purely in terms of $\theta$, which is $$F=\frac{2mMG\sin(\theta)}{d^2+d^2\tan^2(\theta)}$$
Now, given some position (the blue dotted line) $p(t)$, we can set the force equal to $m p''(t)$, which we will have to derive to make it in terms of $\theta$ (as a function of time). I will write stuff out in terms of functions of time to make things more clear. We have $$p(t) = d\tan(\theta(t))$$ so taking the second time derivative yields $$p''(t) = d\sec^2(\theta(t))(\theta''(t) + 2\theta'^2(t)\tan(\theta(t))$$
Since in our setup, theta's acceleration is negative, we have $$\frac{2MG\sin(\theta(t))}{d^2\sec^2(\theta(t))} = -d\sec^2(\theta(t))(\theta''(t) + 2\theta'^2(t)\tan(\theta(t)))$$
Is this differential equation correct? If I set all the constants equal to $1$, the plots Wolfram Alpha outputs make some sort of sense to me, and resemble a harmonic oscillator, but I'm not entirely sure.
Also, is there a simpler way to derive this diffeq given that it is right in the first place? (perhaps, using other formulations of classical mechanics?)
Thanks!


This appears to all be correct. However, I think using $\theta$ as the coordinate has complicated things considerably. More natural is to use $p$ itself as the coordinate, giving $$ m p'' = -\frac{2GMm}{r^2}\frac{p}{r}=- \frac{2GMm p}{(d^2+p^2)^{3/2}} $$ The $3/2$ power in the denominator is annoying, but at least there's neither a first derivative term nor a pile of trig functions. It also makes it easier to spot something important: the RHS is the derivative with respect to $p$ of $2GMm/(d^2+p^2)^{1/2}$. Now, note that $(f(p))' = f'(p)p'$ and $p'p'' = (p'^2/2)'$. If we multiply both sides by $p'$, we can then rewrite the equation as $$ \frac{d}{dt}\left[\frac{mp'^2}{2}-\frac{2GMm}{(d^2+p^2)^{1/2}}\right] = 0. $$ This brings us to the alternative derivation method: the conservation of energy. Indeed, whenever you're working with a conservative force like gravity, the conservation of energy is usually the easier way to get the equations of motion. It has the added benefit that the resulting equation is first-order instead of second-order, though it will have an extra constant for the total energy that you can get from initial conditions.