Simple case [of three-body problem]
The position of the idealized planet in the case shown satisfies the differential equation
∂t, tz[t] -z[t]/(z[t]2 + (1/2 (1 + e Sin[2 π t] ))2 )3/2
where e is the eccentricity of the elliptical orbit of the stars (e = 0.1 in the picture). (Note that the physical situation is unstable: if the planet is perturbed so that there is a difference between its distance to each star, this will tend to increase.) Except when e = 0, the equation has no solution in terms of standard mathematical functions. It can be solved numerically in Mathematica using NDSolve, although a working precision of 40 decimal digits was used to obtain the results shown. Following work by Kirill Sitnikov in 1960 and by Vladimir Alekseev in 1968, it was established that with suitably chosen initial conditions, the equation yields any sequence Floor[t[i + 1] - t[i]] of successive zero-crossing times t[i]. The pictures below show the dependence of z[t] on t and z[0]. As t increases, z[t] typically begins to vary more rapidly with z[0]—reflecting sensitive dependence on initial conditions.