function [upstroke,t,States] = run_single(eNa,eK,gNa,gK,s1,s2)

T = 500; % total simulation time
dt = 0.02; % the time step

if (nargin < 4)
  disp('Not all constants given, will use defalt values:');
  eNa = 50	
  eK = -85
  gNa = 4
  gK = 0.35
  s1 = 50
  s2 = s1;
end

if (nargin == 4)
  disp('No stimulus time given, will use defalt value:');
  s1 = 50
  s2 = s1;
end

if (nargin == 5)
  disp('Running with one stimulus only');
  s2 = s1;
end

stim_time = [s1,s2];

t = 0:dt:T;
States = zeros(length(t),8);

y0 =  [-84,0.5e-7,0,0,1,1,0,1]; % = [V,Ca,x,m,h,j,d,f]
S = zeros(1,1,8);
for i = 1:8
   S(1,1,i) = y0(i);
   States(1,i)= y0(i);
end

V(1) = y0(1);
i_stim = zeros(1,1,8);

for i = 2:length(t)
   i_stim(1,1,1) = 0; 
   for j = 1:length(stim_time)   
     if (t(i) > stim_time(j) && t(i) < stim_time(j)+0.5)
       i_stim(1,1,1) = 80;
     end
   end 
   S = S + dt*ds_dt(t(i), S,[],eNa,eK,gNa,gK)+ dt*i_stim; % Solve the ODEs

   for (j = 1:8)
     States(i,j) = S(1,1,j);
   end
end	   
plot(t,States(:,1));  %xlabel('ms'); ylabel('mV');

d = diff(States(:,1)')./diff(t);
upstroke = max(d);
