function [y,t, derivs]=rinzel

%------------- default values --------------
 
global Cm gNa gK VNa VK s gL VL am Vm aw Vw lambda ton toff Ampl;


Ampl=10;
Cm = 1;
gNa = 100;
VNa = 55;
VK = -72;
s = 1;
gL = 0.3;
gK = 20;
VL = -60;  
am = 0.055;
Vm = -30;
aw = 0.045;
Vw = -47;
lambda = 0.02;


%------------- time period-------------

tspan = [0; 300];          

%--------------- Pulso ----------------

ton=150
toff=250 

%---------   V - y(1)    W - y(2)   X - y(3)  C - y(4) 


%------------- initial values-------------

y0=[  -67.5237 ; 0.1351];
[t,y] = ode45(@f,tspan,y0);
derivs=[];
size(y,1)
for j=1:(size(y,1))
    derivs=[derivs; f(t(j),y(j,:))'];
end
%------------- plot a graph-------------
%  
%hold on;
figure(1)
subplot(2,2,1),
plot(t,y(:,1)),
title(['V']);

subplot(2,2,2),
plot(t,INa(y(:,1),y(:,2)),'m'),
title(['INa']);

subplot(2,2,3),
plot(t,IK(y(:,1),y(:,2)),'c'),
title(['IK']);


%----------------ODE description--------------------

%---------------equations (1), (8) but 1/tau, (14), (17)-------------------

function dydt = f(t,y)
     global ton toff Cm;
     dydt=zeros(2,1);
     dydt(1) =  (1/Cm)*(Iapp(t)-Itot(y(1), y(2)));
     dydt(2) = (W_inf(y(1))-y(2))*tau(y(1)) ; 
 

%----------------equation (10), but 1/tau-----------------------------

function tau1 = tau(z)
     global Vw aw lambda;
     tau1=lambda*(exp(aw*(z-Vw))+exp(-aw*(z-Vw)));
  
function F_inf1 = F_inf(z, a, b)
     F_inf1=1./(1+exp(-2*a.*(z-b)));   

function m_inf1 = m_inf(z)
     global am Vm;
     m_inf1= F_inf(z, am, Vm); 

function W_inf1 = W_inf(z)
     global aw Vw;
     W_inf1= F_inf(z, aw, Vw); 
 
function Itot1 = Itot(a,b)
     Itot1= INa(a,b)+IK(a,b)+IL(a);
     
function INa1 = INa(a,b)
     global gNa VNa;
     INa1= gNa*(m_inf(a)).^3.*(1-b).*(a-VNa);

function IK1 = IK(a,b)
     global gK s VK;
     IK1= gK*(b./s).^4.*(a-VK);

function IL1 = IL(a)
     global gL VL;
     IL1= gL*(a-VL);
     
function Iap_t=Iapp(x)
     global ton toff Ampl;
     if x>ton && x<toff
       Iap_t=Ampl;
     else
       Iap_t=0;
     end