Why is my two body model not working? I have created a simple 1 Dimensional two body particle model in C++.  In the model, particle 1 starts at position (0,0) and particle 2 starts at position (1,0).  The particles are accelerated towards each other at a rate proportional to the distance between them.  
The code looks OK to me, but the results look very bazaar and completely nonphysical.  I suspect the problem may be related to the discrete time steps taken.  I have tried making the rate of acceleration proportional to distance^2 and the results don't look any more realistic.
Here is the code and results:


The above graph is a Position Vs. Time graph.
How can I adjust the code to obtain more reasonable results?
 A: In situations like this, it is a good idea to adjust the time step based on the gradient of the force - because the whole concept of numerical integration is that "things don't change too much from now until the next time step", and that assumption is violated when you move rapidly through a region with fast-changing force.
This has a risky side-effect: if you make the time step "small enough" to cope with the rapid variation, and the variation becomes "infinitely rapid", your algorithm may get "stuck" - you may run into Xeno's paradox of the Achilles and the tortoise.
To avoid getting stuck, some numerical integrations (like Runge-Kutta) are better at taking into account higher order curvature - allowing a bigger step without losing accuracy.
In your case though, the acceleration for the entire next time step is determined by the current position - so if your particles happen to be close to each other at the start of the step, they get thrown far away, and since their attractive force will then be greatly diminished they will have a hard time getting back.
A few criticisms on your code:


*

*You define particle1[4] but only use two components 

*You hard code the time step instead of using a variable like dt

*You use the most basic integration method... please learn about others

*You implicitly define mass = 1 and force constant = 1; consider making those variables (even if you set them to 1)


Regarding the first point, I would use a struct for my particles:
typedef struct{
  double x;
  double v;
} PARTICLE;

and then you can access their properties with
PARTICLE p1, p2;
p1.v = 0.;
p1.x = 0.;
p2.v = 0.;
p2.x = 1.;

Which is immediately more readable... If you want to work in two dimensions, you can either write them explicitly as part of your struct: (double x; double y; double vx; double vy;) or make them arrays.
But really - variable time steps and higher order interpolations would help with the stability of your solution.
A: You're getting particles that get 0 distance from each other, and therefore have infinite acceleration. The computer cannot resolve this, and so it explodes according to how close the step gets to zero. 
I'm not sure what "physical" behaviour would be, given that zero distance point masses are difficult to study. But if you added in a y coordinate and gave it an initial velocity you might get better results (basically an elongated ellipse instead of an oscillating line). 
You could make the initial y velocity very small, so that you approximate the X movement. Or maybe just add a small offset to your division. 
