%
% hhsym0.m
%
% this script computes the lengthy derivatives
% required to obtain the quasiactive model
%
% this requires Matlab's Symbolic Toolbox
%


syms v
v = v - 65;     % stems from use of D&A functionals
am = (v+40)./(10*(1-exp(-(v+40)/10)));
bm = 4*exp(-(v+65)/18);
taum = 1./(am+bm);
minf = am.*taum;

ah = 0.07*exp(-(v+65)/20);
bh = 1./(1+exp(-(v+35)/10));
tauh = 1./(ah+bh);
hinf = ah.*tauh;

an = (v+55)./(100*(1-exp(-(v+55)/10)));
bn = exp(-(v+65)/80)/8;
taun = 1./(an+bn);
ninf = an.*taun;

vr = 0;
dm = diff(minf);
dmr = subs(dm,'v',vr);
dh = diff(hinf);
dhr = subs(dh,'v',vr);
dn = diff(ninf);
dnr = subs(dn,'v',vr);

vr = 0;
bm = subs(minf,'v',vr);
bh = subs(hinf,'v',vr);
bn = subs(ninf,'v',vr);

tmr = subs(taum,'v',vr);
thr = subs(tauh,'v',vr);
tnr = subs(taun,'v',vr);

C = 1;
ENa = 127;
EK = -6;
EL = 3;
EL = 2.8417;

GNa = 120;
GK = 36;
GL = 0.3;

A = [-bm^3*bh*GNa-bn^4*GK-GL  3*bm^2*bh*ENa*GNa  bm^3*ENa*GNa  4*bn^3*EK*GK
     dmr/tmr  -1/tmr 0  0 
     dhr/thr  0  -1/thr 0
     dnr/tnr  0  0  -1/tnr]


eig(A)