ex15ch1.m

these codes are for solving OED with matlab

function ex15ch1
a = 2;
c = 1;
h = input('Specify the step size h (0 < h < 1): ');
if (h <= 0) | (h >= 1)
    error('Must have 0 < h < 1.')
end
% Number of steps to integrate from t = 0 to t = 10.
nsteps = round(10/h);
% Initialize the output arrays.
t = zeros(nsteps,1);
x = zeros(nsteps,1);
y = zeros(nsteps,1);
G = zeros(nsteps,1);
% Initialize for the integration.
tn = 0;
yn = [1; 3];
ynp1 = yn;
% perform the integration.
for i = 1:nsteps
    tnp1 = tn + h;
    ynp1 = yn + h*odes(tn,yn,a,c);
    t(i) = tnp1;
    x(i) = ynp1(1);
    y(i) = ynp1(2);
    G(i) = exp(c*x(i)+a*y(i)) / (x(i)^c * y(i)^a);
    tn = tnp1;
    yn = ynp1;
end
figure
plot(x,y)
xlabel('x(t)')
ylabel('y(t)')
axis([0 4.5 0 3.5])
title('Solution of Volterra''s predator-prey model.')
figure
plot(t,G)
axis([0 10 -Inf Inf])
title('Conserved Quantity G(t,x(t),y(t)).')
%===============================================
function dydt = odes(t,ysol,a,c)
x = ysol(1);
y = ysol(2);
dxdt =  a * (x - x*y);
dydt = -c * (y - x*y);
dydt = [dxdt; dydt];