function three_body_planar ()

  global G m1 m2 m3

  G = 6.67384e-11;

  % Masse Sonne, Erde, Jupiter x 500 [kg]
  m1 = 2e30;
  m2 = 6e24;
  m3 = 500*2e27;

  % AWB Ortsposition [m]
  x01 = [0; 0];
  x02 = [0; 1.5e11];
  x03 = [5*1.5e11; 0];

  % AWB Geschwindigkeiten [m/s]
  v01 = [0; 0];
  tmp = 2*pi*x02(2)/(365*24*3600);
  v02 = [-tmp; 0];
  tmp = 2*pi*x03(1)/(12*365*24*3600);
  v03 = [0; tmp];

  % AWB in Vektoren zusammenfassen
  x0 = [x01(1); x02(1); x03(1); x01(2); x02(2); x03(2); ...
        v01(1); v02(1); v03(1); v01(2); v02(2); v03(2)];

  tend = 20*365*24*3600;

  %'InitialStep' and/or 'MaxStep'
  options = odeset ("RelTol", 1e-9, "AbsTol", 1e-9);
  [t, y] = ode45 (@threebodies, [0, tend], x0', options);

  h = figure;
  hax = axes ();
  hold (hax, "on");
  grid (hax, "on");

  % Bahnkurven
  plot (y(:, 1), y(:, 4), 'r-');
  plot (y(:, 2), y(:, 5), 'g-');
  plot (y(:, 3), y(:, 6), 'b-');

  % Plotten der AWB Massen
  plot (x01(1), x01(2), 'r*');
  plot (x02(1), x02(2), 'g*');
  plot (x03(1), x03(2), 'b*');

  xrange = [min(min(y(:, 1:3))), max(max(y(:, 1:3)))];
  yrange = [min(min(y(:, 4:6))), max(max(y(:, 4:6)))];
  xlim (xrange);
  ylim (yrange);

  %dz = 1/sqrt((xrange(2)-xrange(1))^2+(yrange(2)-yrange(1))^2);
  %daspect ([1, 1, dz]);
  axis equal;

  pause (3);

  for i = 1 : 4 : numel (t)
     tact = tend * (i - 1) / numel(t);
     x1act = interp1 (t, y(:,1), tact);
     y1act = interp1 (t, y(:,4), tact);
     x2act = interp1 (t, y(:,2), tact);
     y2act = interp1 (t, y(:,5), tact);
     x3act = interp1 (t, y(:,3), tact);
     y3act = interp1 (t, y(:,6), tact);
     cla (hax);
     scatter (hax, x1act, y1act, 30, "filled", "o", "markerfacecolor", "r");
     scatter (hax, x2act, y2act, 30, "filled", "o", "markerfacecolor", "g");
     scatter (hax, x3act, y3act, 30, "filled", "o", "markerfacecolor", "b");
     drawnow;
  end

end

% dx = [x1', x2', x3', y1', y2', y3', v1', v2', v3', w1', w2', w3']
% x = [x1, x2, x3, y1, y2, y3, v1, v2, v3, w1, w2, w3]
function dx = threebodies (t, x)

  global G m1 m2 m3

  % Abstand
  dc12 = sqrt ((x(1)-x(2)).^2 + (x(4)-x(5)).^2).^3;
  dc13 = sqrt ((x(1)-x(3)).^2 + (x(4)-x(6)).^2).^3;
  dc23 = sqrt ((x(2)-x(3)).^2 + (x(5)-x(6)).^2).^3;

  % 12 Gleichungen erster Ordnung
  dx = zeros(12, 1);

  % x1' = v1; x2' = v2; x3' = v3
  dx(1) = x(7); dx(2) = x(8); dx(3) = x(9);
  % y1' = w1; y2' = w2; y3' = w3
  dx(4) = x(10); dx(5) = x(11); dx(6) = x(12);

  % v1' = G*m2*(x2-x1)/d12^3 + G*m3(x3-x1)/d13^3
  dx(7) = G * m2 * (x(2)-x(1))./dc12 + G * m3 * (x(3)-x(1))./dc13;
  dx(8) = G * m1 * (x(1)-x(2))./dc12 + G * m3 * (x(3)-x(2))./dc23;
  dx(9) = G * m1 * (x(1)-x(3))./dc13 + G * m2 * (x(2)-x(3))./dc23;

  % w1' = G*m2*(y2-y1)/d12^3 + G*m3(y3-y1)/d13^3
  dx(10) = G * m2 * (x(5)-x(4))./dc12 + G * m3 * (x(6)-x(4))./dc13;
  dx(11) = G * m1 * (x(4)-x(5))./dc12 + G * m3 * (x(6)-x(5))./dc23;
  dx(12) = G * m1 * (x(4)-x(6))./dc13 + G * m2 * (x(5)-x(6))./dc23;

end

