%
%  lect14.m -- ODE test script
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  two ODE choices & IVs with 
%  final time, method (alpha string of .m function file) & options

%Xode = 'Xcubic';  IVs = [1e-2 0];  T = 100;  method = 'ode45';
%options = odeset('reltol',a*1e-3,'abstol',a*1e-6,'refine',1, ...
%		'maxstep',T/10,'stats','on');

Xode = 'XvdPol';  IVs = [           2 0];  T = 400;  method = 'ode15s';
options = odeset('reltol',a*1e-3,'abstol',a*1e-6,'refine',1,'stats','on');

disp(' ')
disp([method ' for ' Xode])
disp(' ')

%  ODE solve with operation count

flops(0);
[t y] = feval(method,Xode,[0 T],IVs,options);
disp(' ')
opcount = flops;  disp(['flops = ' num2str(opcount)])

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  solution plot

Dt = diff(t);  minDt = min(Dt);  maxDt = max(Dt);

disp(' ')
disp(['minDt  = ' num2str(minDt) '  ;  maxDt = ' num2str(maxDt) ...
                                 '  ;  meanDt = ' num2str(mean(Dt))])
disp(['ratio = ' num2str(maxDt/minDt)])

clf;subplot(2,1,1)
plot(t,y(:,1),'kx')
hold on
plot(t,y(:,1),'b')

title('\bf concentrations (blue, red)')
xlabel('\bf t-axis');  ylabel('\bf y-axis');
axis([0 T -2.1 2.1])

subplot(2,1,2)
plot(t(2:end),diff(t),'kx')
title(['\bf timesteps for ' method])
xlabel('\bf t-axis');  ylabel('\bf \Delta t');
axis([0 T 0 1.1*max(diff(t))])

if(Xode == 'Xcubic')
   figure(2);clf
   check = (y(:,2).^2)/2 - (y(:,1).^2)/2 + (y(:,1).^4)/4 ;
   disp(['check = ' num2str(max(check)-min(check))])
   disp(' ')
   plot(t,check,'kx');  hold on
   plot(t,check,'c')
   figure(1)
   subplot(2,1,1)
   plot(t,y(:,2),'kx')
   plot(t,y(:,2),'r')
end