function S=membrane2d
ff=figure;
N=10000;
M=500;
% xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
C_d=ones(M,M);
E_L=ones(M,M)*(-80);
g_L=ones(M,M)*8;
g_Na=ones(M,M)*20;
g_K=ones(M,M)*9;
E_Na=ones(M,M)*60;
g_K_s=ones(M,M)*20;
g_Na_s=ones(M,M)*3;
E_K=ones(M,M)*(-90);
E_K_s=ones(M,M)*(-90);
V=ones(M,M)*(-63.5); 
V_rest=ones(M,M)*(-65); 
I=zeros(M,M);
tau = 200/N; tspan = 0:tau:200;
i=1;
Ind=zeros(M*M*5,1);
Jnd=zeros(M*M*5,1);
Val=zeros(M*M*5,1);
s=1;
disp('Creating sparse matrix, this may take a while...');
drawnow;
for i=2:M-1
    fprintf('%d/100\r',floor(100*i/(M))); drawnow;
    j=1;
    v=-3;Ind(s)=(i-1)*M+j; Jnd(s)=(i-1)*M+j; Val(s)=v; s=s+1;
    Ind(s)=(i-1)*M+j; Jnd(s)=(i-2)*M+j; Val(s)=1; s=s+1;
    Ind(s)=(i-1)*M+j; Jnd(s)=(i)*M+j; Val(s)=1; s=s+1;
    Ind(s)=(i-1)*M+j+1; Jnd(s)=(i-1)*M+j; Val(s)=1; s=s+1;         
    for j=2:M-1
       v=-4; Ind(s)=(i-1)*M+j; Jnd(s)=(i-1)*M+j; Val(s)=v; s=s+1;
       Ind(s)=(i-1)*M+j; Jnd(s)=(i-2)*M+j; Val(s)=1; s=s+1;
       Ind(s)=(i-1)*M+j; Jnd(s)=(i)*M+j; Val(s)=1; s=s+1;
       Ind(s)=(i-1)*M+j-1; Jnd(s)=(i-1)*M+j; Val(s)=1; s=s+1;
       Ind(s)=(i-1)*M+j+1; Jnd(s)=(i-1)*M+j; Val(s)=1; s=s+1;      
    end;
    j=M;
    v=-3;Ind(s)=(i-1)*M+j;Jnd(s)=(i-1)*M+j;Val(s)=v;s=s+1;
    Ind(s)=(i-1)*M+j;Jnd(s)=(i-2)*M+j;Val(s)=1;s=s+1;
    Ind(s)=(i-1)*M+j;Jnd(s)=(i)*M+j;Val(s)=1;s=s+1;
    Ind(s)=(i-1)*M+j-1;Jnd(s)=(i-1)*M+j;Val(s)=1;s=s+1; 
end;
i=1;
for j=1:M
    v=-1;
    if (j>1) v=v-1; end;
    if (j<M) v=v-1; end;
    Ind(s)=(i-1)*M+j;Jnd(s)=(i-1)*M+j;Val(s)=v;s=s+1;
    Ind(s)=(i-1)*M+j;Jnd(s)=(i)*M+j;Val(s)=1;s=s+1;
    if (j>1)
       Ind(s)=(i-1)*M+j-1;Jnd(s)=(i-1)*M+j;Val(s)=1;s=s+1;
    end;
    if (j<M)
       Ind(s)=(i-1)*M+j+1;Jnd(s)=(i-1)*M+j;Val(s)=1;s=s+1;
    end;   
end;    
i=M;
for j=1:M
    v=-1;
    if (j>1) v=v-1; end;
    if (j<M) v=v-1; end;   
    Ind(s)=(i-1)*M+j;Jnd(s)=(i-1)*M+j;Val(s)=v;s=s+1;
    Ind(s)=(i-1)*M+j;Jnd(s)=(i-2)*M+j;Val(s)=1;s=s+1;
    if (j>1)
       Ind(s)=(i-1)*M+j-1;Jnd(s)=(i-1)*M+j;Val(s)=1;s=s+1;
    end;
    if (j<M)
       Ind(s)=(i-1)*M+j+1;Jnd(s)=(i-1)*M+j;Val(s)=1;s=s+1;
    end;   
end;
Ind=Ind(1:s-1);
Jnd=Jnd(1:s-1);
Val=Val(1:s-1);
S=2*sparse(Ind,Jnd,Val);
disp('Done');
n=n_inf(V);
s=s_inf(V);
s1=s1_inf(V);
I=zeros(M,M);
for t=tspan
    I(:)=0;
    if (t>1)&&(t<1.5) I((M/4),(M/4))=500; end; 
    % to ignite any interesting action 
    V = V + tau*((randn(M,M)*5 + I + reshape(S*V(:),M,M) - ...
        g_L.*(V-E_L)-g_Na.*m_inf(V).*(V-E_Na)-g_Na_s.*s1.*(V-E_Na)- ...
        g_K.*n.*(V-E_K)-g_K_s.*s.*(V-E_K_s))./C_d);
    n = n + tau*(n_inf(V) - n);    
    s = s + tau*(s_inf(V) - s)./50; 
    s1 = s1 + tau*(s1_inf(V) - s1)./20;
    i=i+1;
    if (mod(i,10)==0)
       cla; imagesc(V);
       disp(num2str(t));
       colorbar;drawnow;
       % Uncomment to print frames of the animation:
       print(ff, '-r150', '-dpng', ['burst2d', num2str(i),'.png']);
    end;
end;
return;


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

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

function ret = s_inf(v)
   ret=1./(1+exp((-20-v)./5));
   
function ret = s1_inf(v)
   ret=1./(1+exp((-40-v)./5));




