% m29b_planets.m  Orbits of three planets going off to infinity
%
% Here we use one of MATLAB's adaptive ODE solvers to track the orbits
% of three planets attracting each other with pairwise 1/r^2 forces
% initially at rest at the corners of a 3-4-5 right triangle.  The
% system is transiently chaotic: at around t=86, it "self-ionizes".
% Complex arithmetic is used for brevity.
%
% For a discussion of this problem, and an equivalent Chebfun code,
% see www.chebfun.org/examples/ode-non/ThreePlanets.html.

% Set the mood by plotting some "stars":
  function m29()
  fill(20*[-1 1 1 -1 -1],20*[-1 -1 1 1 -1],'k')
  hold on, grid on, axis([1.27*[-6.3 4.7] -3.5 7.5])
  x = 15*rand(250,1)-8; y = 11*rand(250,1)-3.5;
  plot(x,y,'.w','markersize',4)

% Initial conditions:
  x = 0; y = 3.001; z = 4i; u0 = [x y z 0 0 0].';
  xh = plot(real(x),imag(x),'.r','markersize',50);
  yh = plot(real(y),imag(y),'.y','markersize',50);
  zh = plot(real(z),imag(z),'.g','markersize',50);
  tol = 3e-14; 
  opts = odeset('reltol',tol,'abstol',tol);
  axis off, title('t = 0','fontsize',18), drawnow, shg, pause, shg

% Time-stepping:
  dt = .4; tmax = 200;
  tspan = 0:dt:tmax;
  [t,u] = ode113(@m29fun,tspan,u0,opts);
  for nt = 2:length(t);
    tnext = nt*dt;
    x = u(nt,1); y = u(nt,2); z = u(nt,3);
    set(xh,'xdata',real(x),'ydata',imag(x))
    set(yh,'xdata',real(y),'ydata',imag(y))
    set(zh,'xdata',real(z),'ydata',imag(z))
    errest = tol*10^(1+tnext/12);
    title(sprintf('t = %3.0f        estimated error = %5.0e',...
      tnext,errest),'fontsize',18), shg, pause(.06), shg, drawnow
  end

  function v = m29fun(t,u)
  x = u(1); y = u(2); z = u(3);
  yx = (y-x)/abs(y-x)^3;
  zx = (z-x)/abs(z-x)^3;
  zy = (z-y)/abs(z-y)^3;
  v = [u(4:6); yx+zx; -yx+zy; -zx-zy];