Page 220
% first create 12 m-files
%   function y=alphan(v); % (1) alphan.m
%   vrest=-69;
%   y=0.01*(10-(v-vrest))/(exp(0.1*(10-(v-vrest)))-1);
%   function y=betan(v); % (2) betan.m
%   vrest=-69; y=0.125*exp(-(v-vrest)/80);
%   function y=alpham(v); % (3) alpham.m
%   vrest=-69;
%   y=0.1*(25-(v-vrest))/(exp(0.1*(25-(v-vrest)))-1);
%   function y=betam(v); % (4) betam.m
%   vrest=-69; y=4*exp(-(v-vrest)/18);
%   function y=alphah(v); % (5) alphah.m
%   vrest=-69; y=0.07*exp(-0.05*(v-vrest));
%   function y=betah(v); % (6) betah.m
%   vrest=-69; y=1/(exp(0.1*(30-(v-vrest)))+1);
%   function y=pulse(t); % (7) pulse.m
%   if t<.001
%     y=0;
%   elseif t<.002
%     y=-20;
%   else
%     y=0;
%   end
%   function y=rhsV(t,V,n,m,h); % (8) rhsV.m
% Ena=55; Ek=-82; El= -59; gkbar=24.34; gnabar=70.7;
% gl=0.3; cm=0.001;
% y=-(gnabar*m^3*h*(V-Ena)+
%   gkbar*n^4*(V-Ek)+gl*(V-El)+pulse(t))/cm;
%   function y=rhsn(t,V,n,m,h); % (9) rhsn.m
%   y=1000*(alphan(V)*(1-n) - betan(V)*n);
%   function y=rhsm(t,V,n,m,h); % (10) rhsm.m
%   y=1000*(alpham(V)*(1-m) - betam(V)*m);
%   function y=rhsh(t,V,n,m,h); % (11) rhsh.m
%   y=1000*(alphah(V)*(1-h) - betah(V)*h);
%   function Yprime=hhRHS(t,Y); % (12) hhRHS.m
% Yprime=[rhsV(t,Y(1),Y(2),Y(3),Y(4));rhsn(t,Y(1),Y(2),Y(3),Y(4));
%    rhsm(t,Y(1),Y(2),Y(3),Y(4));rhsh(t,Y(1),Y(2),Y(3),Y(4))];
%
vrest=-69;
[t,sol]=ode45('hhRHS',[0 .02], [vrest; 0.315; 0.042; 0.608]);
Vs=sol(:,1);
plot(t,Vs)