function ret = burster()
figure;
N=150000;
C_d=1;
E_L=-80;
I=0;
tau_V=0.152;
g_L=8;
g_Na=20;
g_K=9;
g_K_s=5;
E_Na=60;
E_K=-90;
E_K_s=-90;
V=-65;  
n=n_inf(V);
s=s_inf(V);
VV=zeros(N,1);  II=zeros(N,1); nn=zeros(N,1); ss=zeros(N,1);
tau = 200/N; tspan = 0:tau:200;
i=1;
for t=tspan
    I=5;
    if (t>10)&&(t<12) I=10; end;
    V = V + tau*((I - g_L*(V-E_L)-g_Na*m_inf(V)*(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)/0.152;    
    s = s + tau*(s_inf(V) - s)/20;    
    VV(i)=V;
    nn(i)=n;
    ss(i)=s;
    II(i)=I;
    i=i+1;
end;
hold on;
h=plot(tspan,VV,tspan,-80+II,'linewidth',2);
plot(tspan,ss*100-85,'r','linewidth',2);
axis([0 max(tspan) -90 -5])
ylabel('V');
xlabel('Time (ms)');
legend('Neuron response','Input current','Slow current');
figure
plot3(nn,tspan,VV,'Linewidth',1.5);
xlabel 'n';
ylabel 'time';
zlabel 'V';

figure
h=plot(tspan,ss);
ylabel('Slow current');
xlabel 'time';

figure
plot3(nn,ss,VV,'Linewidth',1.5);
xlabel 'n';
ylabel 'Slow current';
zlabel 'V';

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



