clear
reset(0)
echo on
%This program illustrates the use of an adaptive
%step size algorithm to improve the accuracy of
%integration.  Basically, we want to put in more
%points in the region where the function curvature
%is large, much the same way as was the case in
%adaptive quadrature.  In this case we will use the
%adaptive quadrature routine ode23.m.  We pick
%as our example the predator - prey system that was
%discussed in class.  We shall take the rabbit population
%to be given by y(1) and the fox population to be y(2).
pause
%The populations vary according to the equations:
%f(1)=a*y(1)-b*y(1)*y(2);
%f(2)=(c*y(1)-1)*y(2);
%where a=2, b=1, and c=2 are some constants describing
%the growth and coupling of the populations.  This
%information will be stored in the attached program
%rabbit.m.
pause
%We must set initial conditions, and decide on a tolerance.
%Thus:
options = odeset('RelTol',1e-5);
y0=[1,1];
%We want to integrate this from t=0 to t=15.  Thus:
pause
[t,y]=ode23('rabbit',[0,15],y0,options);
plot(t,y)
set(gca,'FontSize',18)
xlabel('time')
ylabel('fox and rabbit populations')
legend('rabbits','foxes')
pause
%As you can see, the function is rather cyclical with
%first the fox population growing, and then the rabbit
%population exploding after the foxes have died off.  We
%can examine the evolution of the step size as well:
n=length(t);
dt=t(2:n)-t(1:n-1);
hold on
plot(t(1:n-1),dt*40,'w')
set(gca,'FontSize',18)
legend('rabbits','foxes','step size')
pause
%As you can see, there is a tremendous decrease in the step
%size in the regions where either of the functions have a
%high degree of curvature.
hold off
echo off


________________________________


function f=rabbit(t,y)
a=2;
b=1;
c=2;
f=zeros(size(y));
f(1)=a*y(1)-b*y(1)*y(2);
f(2)=(c*y(1)-1)*y(2);