%% Set parameters
v = 2;
c = 1;
ode_fun = @(t,y) hund(t, y, v, c);

y0=[0; 5];
t_span=[0,10];

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

%% Solve ODE by various ode solvers

%tic;[t45,y45]=ode45(ode_fun,t_span,y0,options);time45=toc;
%tic;[t23,y23]=ode23(ode_fun,t_span,y0,options);time23=toc;
%tic;[t113,y113]=ode113(ode_fun,t_span,y0,options);time113=toc;
tic;[t15s,y15s,te15s,ye15s]=ode15s(ode_fun,t_span,y0,options);time15s=toc;
tic;[t23s,y23s,te23s,ye23s]=ode23s(ode_fun,t_span,y0,options);time23s=toc;
tic;[t23t,y23t,te23t,ye23t,ie23t]=ode23t(ode_fun,t_span,y0,options);time23t=toc;
tic;[t23tb,y23tb,te23tb,ye23tb]=ode23tb(ode_fun,t_span,y0,options);time23tb=toc;

%% Plot results

figure;
subplot(2,2,1);plot(y15s(:,1),y15s(:,2),'-*');
xlabel(strcat('steps: ',sprintf('%d',size(t15s)),'  time: ',sprintf('%d',time15s)));
title('ode15s');
ylim([0 y0(2)]);
subplot(2,2,2);plot(y23s(:,1),y23s(:,2),'-*');
xlabel(strcat('steps: ',sprintf('%d',size(t23s)),'  time: ',sprintf('%d',time23s)));
title('ode23s');
ylim([0 y0(2)]);
subplot(2,2,3);plot(y23t(:,1),y23t(:,2),'-*');
xlabel(strcat('steps: ',sprintf('%d',size(t23t)),'  time: ',sprintf('%d',time23t)));
title('ode23t');
ylim([0 y0(2)]);
subplot(2,2,4);plot(y23tb(:,1),y23tb(:,2),'-*');
xlabel(strcat('steps: ',sprintf('%d',size(t23tb)),'  time: ',sprintf('%d',time23tb)));
title('ode23tb');
ylim([0 y0(2)]);


figure;
subplot(2,2,1);plot(t15s, y15s,'-*');
xlabel(strcat('steps: ',sprintf('%d',size(t15s)),'  time: ',sprintf('%d',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('%d',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('%d',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('%d',time23tb)));
title('ode23tb');
%axis([0 1000 -0.1 1]);