%% Set parameters
mu=1/82.45;
ode_fun = @(t,y) dreikoerper(t, y, mu);

y0=[1.2; 0 ; 0; -1.05];
t_span=[0,10];

options=[];
options=odeset('RelTol',1e-6,'AbsTol',1e-12);

%% Solve ODE by various ode solvers

%% ODE45  Solve non-stiff differential equations, medium order method.
tic;[t45,y45]=ode45(ode_fun,t_span,y0,options);time45=toc;

%% ODE23  Solve non-stiff differential equations, low order method.
tic;[t23,y23]=ode23(ode_fun,t_span,y0,options);time23=toc;

%% ODE113 Solve non-stiff differential equations, variable order method.
tic;[t113,y113]=ode113(ode_fun,t_span,y0,options);time113=toc;

%% ODE15S Solve stiff differential equations and DAEs, variable order method.
tic;[t15s,y15s]=ode15s(ode_fun,t_span,y0,options);time15s=toc;

%% ODE23S  Solve stiff differential equations, low order method.
tic;[t23s,y23s]=ode23s(ode_fun,t_span,y0,options);time23s=toc;

%% ODE23T  Solve moderately stiff ODEs and DAEs, trapezoidal rule.
tic;[t23t,y23t]=ode23t(ode_fun,t_span,y0,options);time23t=toc;

%% ODE23TB Solve stiff differential equations, low order method.
tic;[t23tb,y23tb]=ode23tb(ode_fun,t_span,y0,options);time23tb=toc;

%% Plot results

figure;
subplot(3,3,1);plot(y45(:,1),y45(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t45)),'  time: ',sprintf('%d',time45)));
title('ode45');
subplot(3,3,2);plot(y23(:,1),y23(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t23)),'  time: ',sprintf('%d',time23)));
title('ode23');
subplot(3,3,3);plot(y113(:,1),y113(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t113)),'  time: ',sprintf('%d',time113)));
title('ode113');
subplot(3,3,4);plot(y15s(:,1),y15s(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t15s)),'  time: ',sprintf('%d',time15s)));
title('ode15s');
subplot(3,3,5);plot(y23s(:,1),y23s(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t23s)),'  time: ',sprintf('%d',time23s)));
title('ode23s');
subplot(3,3,6);plot(y23t(:,1),y23t(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t23t)),'  time: ',sprintf('%d',time23t)));
title('ode23t');
subplot(3,3,7);plot(y23tb(:,1),y23tb(:,2),'--ok');
xlabel(strcat('steps: ',sprintf('%d',size(t23tb)),'  time: ',sprintf('%d',time23tb)));
title('ode23tb');

figure;
subplot(3,3,1);plot(t45,y45(:,3),'-r',t45,y45(:,4),'--b','Linewidth',2);
title('ode45');
subplot(3,3,2);plot(t23,y23(:,3),'-r',t23,y23(:,4),'--b','Linewidth',2);
title('ode23');
subplot(3,3,3);plot(t113,y113(:,3),'-r',t113,y113(:,4),'--b','Linewidth',2);
title('ode113');
subplot(3,3,4);plot(t15s,y15s(:,3),'-r',t15s,y15s(:,4),'--b','Linewidth',2);
title('ode15s');
subplot(3,3,5);plot(t23s,y23s(:,3),'-r',t23s,y23s(:,4),'--b','Linewidth',2);
title('ode23s');
subplot(3,3,6);plot(t23t,y23t(:,3),'-r',t23t,y23t(:,4),'--b','Linewidth',2);
title('ode23t');
subplot(3,3,7);plot(t23tb,y23tb(:,3),'-r',t23tb,y23tb(:,4),'--b','Linewidth',2);
title('ode23tb');
