%calculate the speed of propagation for the dimensionless HH eq
function HH_speed
global gnabar

 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

dv = 0.001;
V0 = [-12:dv:115];
gnabar = 120;  %nominal value

% first, we estimate the slope of I_ion at the rest state

Ion_current = HH_rhs(V0);

 
%find 3 zeros

ndx = find(Ion_current(1:end-1).*Ion_current(2:end)<=0)
   
 % slope at the rest point
 RM_inv = (Ion_current(ndx(1)+1)-Ion_current(ndx(1)))/dv;
 
 RM = 1/RM_inv
  
%Now simulate the cable equation
 
 L = 2.5; % length of the domain
d = 1; % diffusion coefficient
 
 
N=1000;                   % number of spatial grid points
M = 15;
h = L/N;
dt = h/20;
rat = ones(M,3);
 scal = d*dt/h^2;

X = h*(1:N)';
t0 = L/(16*d);
 sc = 0.2*sqrt(4*pi*d*t0);
kstep = 50;
t=0;
tn = zeros(M,1);
 
% set up the CN matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);
 
lam = 0.3;
  V = 100*sech(X/lam).^2; %initial profile
  
  Vs(:,1) = V;  %save the initial profile
 figure(1)
 plot(X,V,'--','linewidth',2)
 hold on 
   
% uses Crank Nicolson to solve the HH cable equation
 % The algorithm

for n=2:M
for k = 1:kstep
% F = V.*(1-V).*(V-alf);% bistable equation; 
% F = V.*(1-V); %Fisher equation
 F = -RM*HH_rhs(V);
V_t =  scal/2 * ([V(2:N );V(N )] - 2*V + [V(1 );V(1:N-1 )]) +dt*F ;
V  = Am\(V+V_t);
t=t + dt;
 
end
   
 
     figure(1)
 plot(X,V,'linewidth',2)
 Vs(:,n) = V;  % save the current solution
  
 hold on
 %pause(0.05)
tn(n) = t;

end
  
tn(3)-tn(2)

 
figure(1)
xlabel('\xi - dimensionless','fontsize',22)
ylabel('V (mV)','fontsize',20)
 hold off
 
  % here is the computation of speed
TL = t;
 figure(3)
sfc=contour(tn,X,Vs,[50 50],'linewidth',2);
sfcl=length(sfc(1,:));

% to get a good estimate of the speed, this curve needs to be truncated
% some
 
spx = sfc(1, 2:end-40);
spy = sfc(2, 2:end-40);
 
  % now use sfc to estimate the speed
  
  p=polyfit(spx,spy,1) %p(1) is the estimated speed
  spest = p(2)+p(1)*spx;
  
  figure(3)
plot(spx,spy,spx,spest,'--')
 ylabel('\xi', 'fontsize',20)
  xlabel('\tau', 'fontsize',20)
 

function Ion_current = HH_rhs(V)
global gnabar
gkbar = 36;

h = 0.5961;
n = 0.3177;

AM=.1*(25.-V)./(exp(.1*(25.-V))-1.);
BM=4.*exp(-V/18);
m = AM./(AM+BM);

gna=gnabar*m.^3.*h;
 
gk=gkbar*n.^4;

Vna = 115 ;
Vk =  -12 ;
Vl =  10.5988;
Icl=.3*(V-Vl);

Ion_current = -gna.*(V-Vna)-gk.*(V-Vk)-Icl;
