function ret = neurons()
% Sample matlab program with a number of coupled
% neurons (persistent sodium + potassium model)
% Good starting point to a student project.
figure;
N=8; % Number of neurons
% Parameters:
C_d=ones(1,N);
E_L=ones(1,N)*(-80)-((1:N)-1)/4;
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)-2.5*((1:N)-1);
n=n_inf(V,nV);
% Some variables
tau = 0.01; tspan = 0:tau:200;
VV=zeros(size(tspan,1),N);% preallocation, works faster! 
i=1;
for t=tspan
    % The input current is composed of inputs from
    % sorrounding neurons + noise
    % The current decays as I*0.98^t;
    I=I*0.98;
    I=I+circshift((V-circshift(V,[0,1])).*0.01,[0,1]);  % Input from neighboring neurons  
    I=I-circshift((V-circshift(V,[0,-1])).*0.01,[0,-1]);
    % The above part can be done via a clever matrix multiplication 
    % (I=A*V for some A) in order to obtain more complex topology,
    % which is perhaps the best way to do it in Matlab
    I=I+randn(1,N)*3; % some noise
    
    if ((t>18)&(t<32)) I(1,N)=50; end; % N'th neurons receives a 4ms pulse in order 
    % to ignite any interesting action
    
    V = V + tau*((I - 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;    
    VV(i,:)=V;
    i=i+1;
end;
hold on;
for i=1:N % plot all waveforms in 3d
  plot3(tspan*0+i,tspan,VV(:,i),'Linewidth',1);
end;
view(90,80);

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));



