%% Set parameters
k = [0.04, 10^4, 3*10^7];
ode_fun = @(t,y) autokat(t, y, k);

y0=[1; 0 ; 0;];
t_span=[0,1000];

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

%% Solve ODE by various ode solvers

% tic;[t45,y45]=ode45(ode_fun,t_span,y0,options);time45=toc;
% >> size(y45)
% 
% ans =
% 
%      6811789           3
% >> time45 
% 
% time45 =
% 
%    85.3587
% 
% tic;[t23,y23]=ode23(ode_fun,t_span,y0,options);time23=toc;
% 
% >> size(y23)
% 
% ans =
% 
%      2249422           3
% 
% >> time23
% 
% time23 =
% 
%    36.7497
tic;[t113,y113]=ode113(ode_fun,t_span,y0,options);time113=toc;
% 
% >> size(y113)
% 
% ans =
% 
%      3521508           3
% 
% >> time113    
% 
% time113 =
% 
%    95.1879
% 
% 
% tic;[t15s,y15s]=ode15s(ode_fun,t_span,y0,options);time15s=toc;
% tic;[t23s,y23s]=ode23s(ode_fun,t_span,y0,options);time23s=toc;
% tic;[t23t,y23t]=ode23t(ode_fun,t_span,y0,options);time23t=toc;
% tic;[t23tb,y23tb]=ode23tb(ode_fun,t_span,y0,options);time23tb=toc;
% 
% % Plot results
% 
% figure(1); clf;
% subplot(2,2,1);plot(t15s,y15s,'-*');
% xlabel(strcat('steps: ',sprintf('%d',size(t15s)),'  time: ',sprintf('%4.2e',time15s)));
% title('ode15s');
% axis([0 1000 -0.1 1]);
% subplot(2,2,2);plot(t23s,y23s,'-*');
% xlabel(strcat('steps: ',sprintf('%d',size(t23s)),'  time: ',sprintf('%4.2e',time23s)));
% title('ode23s');
% axis([0 1000 -0.1 1]);
% subplot(2,2,3);plot(t23t,y23t,'-*');
% xlabel(strcat('steps: ',sprintf('%d',size(t23t)),'  time: ',sprintf('%4.2e',time23t)));
% title('ode23t');
% axis([0 1000 -0.1 1]);
% subplot(2,2,4);plot(t23tb,y23tb,'-*');
% xlabel(strcat('steps: ',sprintf('%d',size(t23tb)),'  time: ',sprintf('%4.2e',time23tb)));
% title('ode23tb');
% axis([0 1000 -0.1 1]);
