%Exa = 1;

switch Exa
case 1
%% script file to solve a logistic differential equation
x0 = 0.2;  % initial value
t0 = 0;    % initial time
Tend = 10; % time interval for which solution has to be computed
n = 50;    % number of intermediate times

%% define the function describing the right hand side
%function y = logF(t,x)
%  y = x-x.*x;
%endfunction

logF = @(t,x)x*(1-x)

% generate vector with the times at which solution is desired
Tmin = t0 + linspace(0,Tend,n); 
vopt = odeset('RelTol', 1e-3, 'AbsTol', 1e-3);
[T,X] = ode45(logF,Tmin,x0,vopt);

plot(T,X)
axis([0,Tend,0,1]); grid on
[T2,X2] = ode45(logF,Tmin,2*x0,vopt);
[T3,X3] = ode45(logF,Tmin,0.2*x0,vopt);

plot(T,X,T2,X2,T3,X3)
xlabel('time t'); ylabel('population'); grid on

case 2
%%%%%%%%%%% Volterra Lotka %%%%%%%%%%%%%%
%function res = VolterraLotka(t,x)
%  c1 = 1; c2 = 2; c3 = 1; c4 = 1;
%  res = [(c1-c2*x(2))*x(1); (c3*x(1)-c4)*x(2)];
%endfunction
c1 = 1; c2 = 2; c3 = 1; c4 = 1;
VolterraLotka = @(t,x)[(c1-c2*x(2))*x(1); (c3*x(1)-c4)*x(2)];

%%%%%% vector field %%%%%%%%%%%%%
x = 0:0.2:2.6;  % define the x values to be examined
y = 0:0.2:2.0;  % define the y values to be examined

n = length(x); m = length(y);
Vx = zeros(n,m); Vy = Vx; % create zero vectors for the vector field

for i = 1:n
  for j = 1:m
    v = VolterraLotka(0,[x(i),y(j)]); % compute the vector
    Vx(i,j) = v(1); Vy(i,j) = v(2);
  end%for
end%for

figure(1);
quiver(x,y,Vx',Vy',3);
axis([min(x),max(x),min(y),max(y)]);
grid on; xlabel('prey'); ylabel('predator');

case 3
%%%%%%%%%% solutions of ODE system
t  = linspace(0,15,100);
[t,XY] = ode45(VolterraLotka,t,[2;1],t);

figure(2);
plot(t,XY)
xlabel('time'); ylabel('predator and prey')
legend('prey','predator'); axis([0,15,0,3]); grid on

figure(3);
plot(XY(:,1),XY(:,2),'r');
axis([min(x),max(x),min(y),max(y)]);
hold on
quiver(x,y,Vx',Vy',2);
grid on; xlabel('prey'); ylabel('predator'); hold off

case 4  % Octave only
%%%%%%%%% speed of compiled code %%%%%%%%%%%
%vopt = odeset('RelTol', 1e-13, 'AbsTol', 1e-13,'NormControl','on');
vopt = odeset('RelTol', 1e-10, 'AbsTol', 1e-10,'NormControl','on');
%tic()
t0 = cputime();
[t,XY] = ode45(VolterraLotka,[0,50],[2;1],vopt);
TimeScript = cputime()-t0
%TimeScript = toc() 

%tic()
t0 = cputime();
%[t,XY] = ode45(@(t,x)VolterraLotkaC(x,t),[0,50],[2,1],vopt);
[t,XY] = ode45(@VolterraLotkaC,[0,50],[2;1],vopt);
TimeCPP = cputime()-t0
%TimeCPP = toc()
ratio = TimeScript/TimeCPP

case 5
%%%%% periodicity
x0 = 2; y0 = 1;  % initial values
[t,XY] = ode45(VolterraLotka,[0,7],[x0;y0],vopt);
% examine only times between 6 and 7
ind = find((t>=6)&(t<=7));
t = t(ind); x = XY(ind,1);
plot(t,x-x0);  % visual test for zero
grid on                                
s = sign((XY(ind,1)-x0).*(XY(ind-1,1)-x0)); % detect sign changes
pos = find(s<0);                            % position of sign change

% use linear interpolation to determine the partial time step
dt = (t(pos)-t(pos-1))*(x(pos-1)-x0)/(x(pos-1)-x(pos))
T = t(pos-1)+dt                                    % estimate the period
[tt,XYt] = ode45(@(t,x)VolterraLotkaC(x,t),[T-dt,T],XY(ind(pos)-1,:));
% compute the value at T
difference = XYt(end,:)-[x0,y0]  % difference with initial values
end%switch
