%Problem: The script does not run
%
% zooplankton-algae interaction
% The effect of fish on the
%#Ref: Scheffer, M. (1991) "Fish and nutrients interplay determines algal biomass: A minimal model" Oikos 62:271-282.
%
%cons=g*Z*(A/(A+h))
%A' = A*r*(1-A/K) - cons + i*(K-A)
%Z' = e*cons - m*Z - F*Z^2/(Z^2+hz^2)%+sigma*randn(1);
%
function [bd1] = zaf
mod=get_model;
% coco settings
prob = coco_prob();
prob = coco_set(prob, 'ode', 'vectorized','on');
prob = coco_set(prob,'EP', 'NSA',true,'BTP',true);
prob = coco_set(prob,'cont','h0',0.0003,'h_min',0.0003,'h_max',0.03,'h_fac_min',0.5,'h_fac_max',2,'PtMX',400,'NPR',10,'NSV',[]);
prob = coco_set(prob,'corr','ItMX',10,'SubItMX',7,'TOL',1e-006,'ResTOL',1e-006,'NTST',50,'NCOL',5);

%run coco 1D try to follow initial orbit
%initial conditions after long stabilization
A=0.800104891213651;
Z=2.75227633950124;
%run one cycle (approximate period of 33 timesteps)
[ts,Ys]=ode45(mod.odehandle,[0 33],[A,Z]);

%initial condition to po
prob = ode_isol2po(prob, '',mod.coco_ode, mod.coco_jac, mod.coco_jacp,ts,Ys, mod.parnames, mod.p0);
bd1 = coco(prob,'1', [], 1,  'F', [0 0.3]);


function g_l_res=get_model
%parameters:
e  = 0.6;           %,-
F  = 0;             %,g m-3 d-1
g  = 0.4;           %,d-1
h  = 0.6;           %,g m-3
hz = 1;             %,g m-3
i  = 0.04;          %,d-1
K  = 10;            %,g m-3
m  = 0.15;          %,d-1
r  = 0.5;           %,d-1
g_l_res.odehandle=@(t,g_X1)curr_ode2(t,g_X1,e,F,g,h,hz,i,K,m,r);
g_l_res.coco_ode=@(g_X1,g_par)curr_ode2(0,g_X1,g_par(1,:),g_par(2,:),g_par(3,:),g_par(4,:),g_par(5,:),g_par(6,:),g_par(7,:),g_par(8,:),g_par(9,:));
g_l_res.coco_jac=@(g_X1,g_par)reshape(i_vcat(-g_par(6,:)-g_par(9,:).*(g_X1(1,:)./g_par(7,:)-1.0)-(g_X1(2,:).*g_par(3,:))./(g_X1(1,:)+g_par(4,:))-(g_X1(1,:).*g_par(9,:))./g_par(7,:)+g_X1(1,:).*g_X1(2,:).*g_par(3,:).*1.0./(g_X1(1,:)+g_par(4,:)).^2,(g_X1(2,:).*g_par(1,:).*g_par(3,:))./(g_X1(1,:)+g_par(4,:))-g_X1(1,:).*g_X1(2,:).*g_par(1,:).*g_par(3,:).*1.0./(g_X1(1,:)+g_par(4,:)).^2,-(g_X1(1,:).*g_par(3,:))./(g_X1(1,:)+g_par(4,:)),-g_par(8,:)-(g_par(2,:).*g_X1(2,:).*2.0)./(g_X1(2,:).^2+g_par(5,:).^2)+g_par(2,:).*g_X1(2,:).^3.*1.0./(g_X1(2,:).^2+g_par(5,:).^2).^2.*2.0+(g_X1(1,:).*g_par(1,:).*g_par(3,:))./(g_X1(1,:)+g_par(4,:))),[2,2,max(size(g_X1,2),size(g_par,2))]);
g_l_res.coco_jacp=@(g_X1,g_par)reshape(i_vcat(0.0,(g_X1(1,:).*g_X1(2,:).*g_par(3,:))./(g_X1(1,:)+g_par(4,:)),0.0,-g_X1(2,:).^2./(g_X1(2,:).^2+g_par(5,:).^2),-(g_X1(1,:).*g_X1(2,:))./(g_X1(1,:)+g_par(4,:)),(g_X1(1,:).*g_X1(2,:).*g_par(1,:))./(g_X1(1,:)+g_par(4,:)),g_X1(1,:).*g_X1(2,:).*g_par(3,:).*1.0./(g_X1(1,:)+g_par(4,:)).^2,-g_X1(1,:).*g_X1(2,:).*g_par(1,:).*g_par(3,:).*1.0./(g_X1(1,:)+g_par(4,:)).^2,0.0,g_par(2,:).*g_X1(2,:).^2.*g_par(5,:).*1.0./(g_X1(2,:).^2+g_par(5,:).^2).^2.*2.0,-g_X1(1,:)+g_par(7,:),0.0,g_par(6,:)+g_X1(1,:).^2.*1.0./g_par(7,:).^2.*g_par(9,:),0.0,0.0,-g_X1(2,:),-g_X1(1,:).*(g_X1(1,:)./g_par(7,:)-1.0),0.0),[2,9,max(size(g_X1,2),size(g_par,2))]);
g_l_res.p0=[e;F;g;h;hz;i;K;m;r];

g_l_res.parnames={'e','F','g','h','hz','i','K','m','r'};

%function created by GRIND
function g_X2=curr_ode2(t,g_X1,e,F,g,h,hz,i,K,m,r)
cons=g.*g_X1(2,:).*(g_X1(1,:)./(g_X1(1,:)+h));
g_X2(1,:)=g_X1(1,:).*r.*(1-g_X1(1,:)./K)-cons+i.*(K-g_X1(1,:));
g_X2(2,:)=e.*cons-m.*g_X1(2,:)-F.*g_X1(2,:).^2./(g_X1(2,:).^2+hz.^2);

function m=i_vcat(varargin)
len=cellfun('size',varargin,2);
len1=len==1;
if any(len1)&&any(~len1)
    flen1=find(len1);
    maxlen=max(len);
    for i=1:length(flen1)
        varargin{flen1(i)}=varargin{flen1(i)}(ones(1,maxlen));
    end;
end;
m=vertcat(varargin{:});
