I am attempting to create an n-dimensional physics simulation in C++ without libraries. I have made collisions between the n-dimensional spheres work and am currently working with the gravitational attraction between objects. As it is, whenever two objects make a close approach, energy seems to not remain constant as the two objects often depart faster than they approached. I have now been attempting to negate this by solving the two-body problem for my use case. Problem is, with my being in Calculus 1 right now, I don't know how to some certain parts of the problem. As well as this, the problem does not seem to be too well documented online. I have been directed from answers saying things like literally "here's a vague answer for a vague question" to literal books that would take an age to parse though, of which I don't want to spend a few months reading a book to simply stop my objects from exploding over time. I feel like the internet is lacking some middle ground between the two. I checked multiple wikipedia pages, and was sort of lead to where I am about to show, and then was told something along the lines of "from here on, you do a simple so-and-so to derive the solution you're looking for." I am coming here to the mathematics stack exchange instead of the physics stack exchange because I believe the part I got stuck on was a mathematics problem. As well as this, I don't want to use recursive approximations for the real values, as I feel they are both computationally expensive and don't get rid of the problem, as if the two objects' masses were increased or their radii decreased (perhaps by multiple orders of magnitude, but my point still stands), then the problem would come back. I would also like to be able to use larger time steps and have similar precision to that of smaller time steps in my simulation. I feel like using calculus would solve my issue, specifically using integrals in some way, as calculus deals with the instantaneous, of which I would like calculated between my time steps to make them more precise.
Problem: Given two objects of mass $m_1$ and $m_2$ going at initial velocities $v_{1i}$ and $v_{2i}$ at positions $x_{1i}$ and $x_{2i}$ respectively, where v and x are n-dimensional vectors, find the final positions and velocities ($x_{1f}$, $x_{2f}$, $v_{1f}$, and $v_{2f}$) of the two objects after a given time $t$.
This mathematics stack exchange answer and the wikipedia page for the 2-body problem suggested that the problem be split up into two parts, being the acceleration and the constant movement of the center of mass of the two objects, so that is what I did. Now that the center of mass is not a problem anymore for the time being, I am left with the following equation:
$${\mu}a={\mu}\frac{d^{2}r}{dt^{2}}=F\left(r\right)$$
Where F(r) is the amount of force at a given distance away, expressed by the following equation:
$$F\left(r\right)=\frac{Gm_{1}m_{2}}{r^{2}}$$
Where $\mu$ is a constant $\frac{m_{1}m_{2}}{m_{1}+m_{2}}$.
I am now attempting to solve for $x\left(x_{i},v_{i},t\right)$ and $v\left(x_{i},v_{i},t\right)$, where $x_i$ and $v_i$ are the relative position and velocity of one object to the center of mass. I believe that this also means that the other object either has equal and opposite values, or different values scaled by some constant that has to do with the mass, perhaps being ${\mu}$.
I have simplified this into the following:
$$a=\frac{d^{2}r}{dt^{2}}=\frac{k}{r^{2}}$$
where k is a constant equal to $\frac{Gm_{1}^{2}m_{2}^{2}}{m_{1}+m_{2}}$.
I attempted rearranging the equation by multiplying both sides by $dt$ twice and $r^2$:
$$\frac{d^{2}r}{dt^{2}}\cdot dt\cdot dt\cdot r^{2}=\frac{k}{r^{2}}\cdot dt\cdot dt\cdot r^{2}$$
$$r^{2}d^{2}r=kdtdt$$
I then tried expanding the left side of the equation. If $d^{2}r=drdr$ similar to how $dt^{2}=dtdt$, then perhaps I could take the integral of both sides and hopefully $v_i$ would pop out as the constant of the first integral and $x_i$ would pop out as the constant in the second integral. I have tried that already, evaluating them separately and checking my work with wolfram alpha (as we haven't gotten to integrals yet in my class), and found that the $x_i$s that popped out cancel each other out, which makes absolutely zero sense in the context of the problem. By this, I conclude that $dx^2$ is different from $d^2x$, and thus should be treated differently.
$$r^{2}\cdot d^{2}r=kdtdt$$
I was hoping that $x_{i}$ and $v_{i}$ would pop out of the double indefinite integration with respect to time as the constants that I need to plug into the original equation.
I am now left with the following:
$$d^{2}r=kdtdt$$
I am confused on how the number of the second derivative is different from the denominator. And does this distinction apply with the first derivative as well? How can I proceed into finding $x\left(x_{i},v_{i},t\right)$ and $v\left(x_{i},v_{i},t\right)$? Is this possible? Is my thought process of $x_i$ and $v_i$ popping out of the integrals accurate? If I took the integrals, would they be with respect to time? How is the integral of $d^2r$ different from that of $dt$ and $dtdt$ and the second integral of $dtdt$?