%  Kurt E. Anderson
%
%  December 23, 2014
%
%  DESCRIPTION:	solves phfn_obe.m with dde23
%
%  OUT:		vector t:  tinit:tinterval:tfinal
%			vector y:  initial values
%			vector p:  parameter values

function y = phode_obe()		% function call

tic;

patchnum = 101;	        % patch number

% PARAMETERS

p = [1; 5; 12; 30; 0.25];
%b       = p(1);        % induction half-maximum
%theta   = p(2);        % strength of threshold
%D0      = p(3);        % baseline dispersal rate of herbivores
%chi     = p(4);        % sensitivity of herbivore dispersal rate to induction level
%m       = p(5);        % mortality
tauh     = 2.5;

% VARIABLES

initvalues(1:patchnum) = 1; % y(1:patchnum), initial values for induction level
initvalues(patchnum+1:2*patchnum) = 1+.1*sin((0:1:patchnum-1).*5*2*pi/(patchnum-1)); % y(patchnum+1:2*patchnum), initial values for herbivore density
tstart = 0;

% TIMES
tfinal = 101;
times = tstart:0.01:tfinal;			% tinit:tinterval:tfinal

% ODE MODEL
options = ddeset('RelTol',1e-3,'AbsTol',1e-6,'OutputFcn',@odeprog,'Events',@odeabort);
sol = dde23('phfn_obe',tauh,initvalues,times,options,p);
%sol = ode23s('phfn_obe',times,initvalues,options,p); % required if time delay is zero.

t = 0:0.01:sol.x(end);
y = deval(sol,t);
y = y';
% Find good location for 2-D output
[~,I] = max(y(1,patchnum+1:2*patchnum));

% OUTPUT

figure(1)               % 3-D output of transient dynamics
subplot(2,1,1)
mesh(0:1:(patchnum-1),t,y(:,1:patchnum));
title('Induction Level')
xlabel('');
ylabel('Time');
zlabel('Level');
axis([0 patchnum-1 0 tfinal-1])

subplot(2,1,2)
mesh(0:1:(patchnum-1),t,y(:,patchnum+1:2*patchnum));
title('Herbivore Density')
xlabel('Location');
ylabel('Time');
zlabel('Density');
axis([0 patchnum-1 0 tfinal-1])

figure(2)               % 2-D output of single location
subplot(2,1,1)
plot(t,y(:,I));
xlabel('Time');
ylabel('Induction Level');

subplot(2,1,2)
plot(t,y(:,patchnum+I));
xlabel('Time');
ylabel('Herbivore Density');

figure(3)               % 2-D output of single time
subplot(2,1,1)
plot(1:patchnum,y(1,1:patchnum));
xlabel('Location');
ylabel('Induction Level');

subplot(2,1,2)
plot(1:patchnum,y(1,patchnum+1:2*patchnum));
xlabel('Location');
ylabel('Herbivore Density');

figure(4)               % 2-D output of single time
subplot(2,1,1)
plot(1:patchnum,y(end,1:patchnum));
xlabel('Location');
ylabel('Induction Level');

subplot(2,1,2)
plot(1:patchnum,y(end,patchnum+1:2*patchnum));
xlabel('Location');
ylabel('Herbivore Density');

y = sol; % output solution to workspace

toc;
end