function ret = cable()
% Sample matlab program simulating
% a piece of an axon via neuronal cable equation
% with I_Na-I_K model of a single membrane compartment.
ff=figure;
N=512; % Number of spatial compartments (spatial resolution)
% Parameters:
C_d=ones(1,N);
E_L=ones(1,N)*(-80);
g_L=ones(1,N)*8;
g_Na=ones(1,N)*20;
g_K=ones(1,N)*10;
E_Na=ones(1,N)*60;
E_K=ones(1,N)*(-90);
V=ones(1,N)*(-65); 
V_rest=ones(1,N)*(-65); 
I=zeros(1,N);
tau_V=ones(1,N);
nV=ones(1,N)*(-25);
n=n_inf(V,nV);
tau = 0.01; tspan = 0:tau:60;
I=zeros(1,N);

% Matrix S is used to simulate the difussion, that it the second order
% derivative with respect to spatial variable
% d^2f/dx^2=(-2f(x)+f(x+1)+f(x-1))/(delta x)^2
S=diag((-2)*ones(1,N),0)+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
S(1,N)=1;
S(N,1)=1;
S=S*10;
S=sparse(S); % making a sparse matrix
% eventually (S*V')' is the difussion factor
i=1;
for t=tspan
    if (t>0) I(1,N/4+N/2)=50; end; % to ignite any interesting action 
    if (t>30) I(1,N/4+N/2)=0; end;
     
    V = V + tau*((I+(S*V')' - g_L.*(V-E_L)-g_Na.*m_inf(V).*(V-E_Na)-g_K.*n.*(V-E_K))./C_d);
    n = n + tau*(n_inf(V,nV) - n)./tau_V;    
    if (mod(i,10)==0)
       plot(1:N,V,1:N,n*50-100,1:N,I/4-120,'Linewidth',2);
       axis([1 N -120 10]);
       legend('V','n','I');
       drawnow;
       % Uncomment to print frames of the animation:
       % print(ff, '-r150', '-dpng', ['cable', num2str(i),'.png']);
    end;
    i=i+1;
    %disp(t);
end;


function ret = m_inf(v)
   ret=1./(1+exp((-20-v)./15));

function ret = n_inf(v, nV)
   ret=1./(1+exp((nV-v)./5));



