%
% solve the 3 cell hhred net
%
% usage    tricirc(g)
%
% where, e.g., g = [4 4 8]
%

function tricirc(g)

n0 = an(0)/(an(0)+bn(0));
s0 = sinf(0);
y0 = [0 n0 s0 0 n0 s0 0 n0 s0];
[t,y] = ode23(@triwork,[0 40],y0,[],g);
plot(t,y(:,[1 4 7]))
return

%

function val = triwork(t,y,p)
val = zeros(9,1);

val(1,1) = -36*y(2)^4*(y(1)+12)...
           -120*minf(y(1))^3*(4/5-y(2))*(y(1)-115)-3/10*(y(1)-10.613)...
           +p(3)*y(9)*(y(1)+10) + 20*exp(-(10/3*t-10/3)^2);

val(2,1) = an(y(1))*(1-y(2))-bn(y(1))*y(2);

val(3,1) = (sinf(y(1)) - y(3))/taus(y(1));

val(4,1) = -36*y(5)^4*(y(4)+12)...
           -120*minf(y(4))^3*(4/5-y(5))*(y(4)-115)-3/10*(y(4)-10.613)...
           +p(1)*y(3)*(y(4)+10);

val(5,1) = an(y(4))*(1-y(5))-bn(y(4))*y(5);

val(6,1) = (sinf(y(4))-y(6))/taus(y(4));

val(7,1) = -36*y(8)^4*(y(7)+12)...
           -120*minf(y(7))^3*(4/5-y(8))*(y(7)-115)-3/10*(y(7)-10.613)...
           +p(2)*y(6)*(y(7)+10);

val(8,1) = an(y(7))*(1-y(8))-bn(y(7))*y(8);

val(9,1) = (sinf(y(7)) - y(9))/taus(y(7));

return

function val = an(v)
val = (1/10-v/100)/(exp(1-v/10)-1);

function val = bn(v)
val = 1/8*exp(-v/80);

function val = minf(v)
val =  1./(1+40*exp(-v/18).*(exp(5/2-v/10)-1)./(25-v));

function val = sinf(v)
val = 1/(1+exp(6-v/5));

function val = taus(v)
val = 18;