ms = 1.9885e30;
mz = 5.97237e24;
G  = 6.6743015e-11;
x0 = 1.47098074e11; 
v0 = 3.0287e4;
g  = G*ms;
solver_type = 'Euler'; %{'Euler' 'default' 'RK2' 'RK4' 'CrankNicolson' 'ImplicitEuler'}

hh = 1/24;
h  = 60*60*24*hh;       % jedna hodina
n  = floor(365.25/hh);  % celkovy pocet kroku simulace
X0 = [x0; 0; 0; v0];
X  = zeros(4,n+1);  X(:,1) = X0;
F  = @(t,Y) [Y(3:4,1); -g*Y(1:2,1)*norm(Y(1:2,1))^(-3)];
% Resime soustavu pro Y(t) = [x;v](t)
% X'(t) = F(t,X(t)) = [v(t); a(t)] kde a(t) = -G*mz*x(t)*norm(x(t))^-3
t = 0;
for i = 1:n+1
  switch lower(solver_type)
    case 'rk2'
      K1 = F(t    ,X(:,i)         );
      K2 = F(t+h/2,X(:,i)+(h/2)*K1);
      X(:,i+1) = X(:,i) + h * K2; 
    case 'rk4'
      K1 = F(t    ,X(:,i)         );
      K2 = F(t+h/2,X(:,i)+(h/2)*K1);
      K3 = F(t+h/2,X(:,i)+(h/2)*K2);
      K4 = F(t+h  ,X(:,i)+ h   *K3);
      X(:,i+1) = X(:,i) + h * (K1+2*K2+2*K3+K4)/6; 
    case 'cranknicolson'
      K1 = F(t    ,X(:,i)         );
      K2 = F(t+h  ,X(:,i)+ h   *K1);
      X(:,i+1) = X(:,i) + h * (K1+K2)/2;
    case 'impliciteuler'
      X(:,i+1) = X(:,i);
      for ii = 1:10
        K = F(t+h,X(:,i+1));
        X(:,i+1) = X(:,i) + h * K;
      end
    otherwise % 'euler'
      K = F(t,X(:,i));
      X(:,i+1) = X(:,i) + h * K; 
  end
  t = t + h;
end
figure(1); hold on;
plot( 0,0,'xr');
plot(x0,0,'xb');
line(X(1,:),X(2,:));
axis equal;