This problem can be treated as an elliptical Kepler orbit. But for a Kepler orbit it is assumed that one mass is much more massive than the other: $m_1>>m_2$, which means that $m_1$ basically sits still at their center of mass.
But both masses will move symmetrically around their center of mass (if one moves closer the other will move closer as well, inversely proportional to their masses), which allows you to write their attractive force as a function of their distance towards their center of mass ($r_{COM}$).
$$
F_1=\frac{Gm_1m_2}{\left(\frac{m_1+m_2}{m_2}r_{COM}\right)^2}
$$
This only scales the force and is equal to an object orbiting a much more massive object (located at the previous center of mass), like a Kepler orbit, but with different masses.
However Kepler orbits do not have an explicit solution for the position as a function of time, these are often calculated numerically. But you can calculate it explicitly the other way around, time as a function of position and the exact trajectory an object will take can also be found explicitly.
Edit:
You can also approximate an orbit with a Fourier series which has the advantage that it goes on infinitely, but will contain small errors. I did some testing with this and got the following results:
$$\theta(t)=\bar{\omega}t+\sum_{n=1}^\infty{A(n)\sin\left(n\bar{\omega}t\right)}$$
for $e=0.5$ (orbital eccentricity) $A(n)\approx 0.9757633n^{-1.944954}$.
$$
v(t)=\sqrt{\frac{\mu}{a(1-e^2)}(1+e(2\cos{\theta(t)}+e))}
$$
When choosing $\frac{\mu}{a}=1$, using the approximation and a limited sum of $n=50$, the graph looks like this:
