%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This page contains MatCont input and MatLab commands to accompany a tutorial
on MatCont. Please look at the instructions file on my website for more 
information. Feel free to use, and adjust. This file page is plain ASCII so
copy and paste is possible.
Last edit: H.G.E. Meijer, 17-2-14
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% JansenRitslow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Name 	    JR_slow
Coordinates y0,y4,y1,y5,y2,y6,y3,y7
Parameters  A,a,Bs,bs,Bf,bf,e0,r,v0,theta,CC,II
Equations
C1=CC;
C2=0.8*CC;
C3=0.25*CC;
C4=0.25*CC;
S4=2*e0/(1+exp(r*(v0-(y1-theta*y2-(1-theta)*y3))));
S5=2*e0/(1+exp(r*(v0-(C1*y0))));
S6=2*e0/(1+exp(r*(v0-(C3*y0))));
S7=S6;
y0'=y4;
y4'=A*a*S4-2*a*y4-a*a*y0;
y1'=y5;
y5'=A*a*(II+C2*S5)-2*a*y5-a*a*y1;
y2'=y6;
y6'=Bf*bf*C4*S6-2*bf*y6-bf*bf*y2;
y3'=y7;
y7'=Bs*bs*C4*S7-2*bs*y7-bs*bs*y3;

%Plot profile
ii=10;jj=1;ndim=8;ntst=40;ncol=4;
mesh=sort([reshape(repmat(f(1:ntst,ii),1,ncol)+...
     repmat((0:ncol-1)/ncol,ntst,1).*repmat(diff(f(1:ntst+1,ii)),1,ncol),ntst*ncol,1); 1]);
plot(mesh*x(end-1,ii),x(jj:ndim:end-2,ii))

%SWD commands
load Systems\JR_slow\diagram\lc1.mat
out=x(3:8:end-2,:)-(x(5:8:end-2,:)+x(7:8:end-2,:))/2;
extrema=nan(6,size(out,2));
for i=1:size(out,2)
diff_tr=diff(out(1:end-1,i));
diff_tr_shift=diff(out(2:end,i));
temp=[out(find([ (diff_tr>0).*(diff_tr_shift<0) ]),i);...
	out(find([ (diff_tr<0).*(diff_tr_shift>0) ]),i)];
extrema(1:size(temp,1),i)=temp;
end
figure(1);hold on;plot(x(end,:),sort(extrema));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% LORENZ-84
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Name 	    Lorenz84
Coordinates x,y,z
Parameters  a,b,F,G
Equations
x'=-y^2-z^2-a*x+a*F
y'=x*y-b*x*z-y+G
z'=b*x*y+x*z-z

Plotting commands:
%%%
load Systems/diagram/H_LC(1).mat
ntst=20;ii=4;jj=40;ndim=3;
figure;plot(f(1:ntst+1,ii)*x(end-1,ii),x(jj:ndim:end-2,ii))
%%%
ntst=20;ii=4;jj=40;ndim=3;ncol=4
plot(f(1:ntst+1,ii)*x(end-1,ii),x(jj:ndim*ncol:end-2,ii))
%%%
ii=240;jj=1;ndim=3;ntst=20;ncol=4;
mesh=sort(reshape(repmat(f(1:ntst,ii),1,ncol)+...
     repmat((0:ncol-1)/ncol,ntst,1).*repmat(diff(f(1:ntst+1,ii)),1,ncol),ntst*ncol,1); 1]);
plot(mesh*x(end-1,ii),x(jj:ndim:end-2,ii))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SEIR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Name 		SEIR
Coordinates ss,ee,ii,x,y
Parameters  alpha,beta0,gamma,delta,mu,omega
Equations
beta=beta0*(1+delta*x)
ss'=mu*exp(-ss)-mu-beta*exp(ii)
ee'=beta*exp(ss+ii-ee)-(mu+alpha)
ii'=alpha*exp(ee-ii)-(mu+gamma)
x'=-omega*y+x*(1-x^2-y^2)
y'=omega*x+y*(1-x^2-y^2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Hodgkin-Huxleuy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Name	    HodgkinHuxley
Coordinates V,m,h,n
Parameters  Istim
Equations
gna=120
gk=36
gleak=.003
Ena=50
Ek=-90
Eleak=-70
am=0.1*(V+35)/(1-exp(-(V+35)/10))
bm=4*exp(-(V+60)/18)
ah=0.07*exp(-(V+60)/20)
bh=1/(1+exp(-(V+30)/10))
an=0.01*(V+50)/(1-exp(-(V+50)/10))
bn=0.125*exp(-(V+60)/80)
minf=am/(am+bm)
hinf=ah/(ah+bh)
ninf=an/(an+bn)
taum=1/(am+bm)
tauh=1/(ah+bh)
taun=1/(an+bn)
V'=-(gna*m^3*h*(V-Ena)+gk*n^4*(V-Ek)+gleak*(V-Eleak))+Istim
m'=(minf-m)/taum
h'=(hinf-h)/tauh
n'=(ninf-n)/taun


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FHN in travelling coordinates
%%% for waves and pulses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Name        FitzhughNagumo_wave
Coordinates U,W,V
Parameters  aa,bb,cc
Equations
U'=W
W'=cc*W+U*(U-aa)*(U-1.0)+V
V'=bb*U/cc