function [output] = DiffEquations(time,init)

global vNa vK vL vH ...
    gNa gK gL gCa gSK gNap gKNa gH gA ...
    tauNbar tauHPbar tauE tauH Ca_ex RTF ...
    thetaM sigmaM thetaN sigmaN thetaS sigmaS ...
    thetaMP sigmaMP thetaHP sigmaHP thetaRF sigmaRF thetaRS sigmaRS ...
    kr prf taurs tauhp f ...
    thetaA sigmaA thetaE sigmaE ...
    C iapp eps kca ks pmax ec nh alphaNa rpump kp neq ...
    gT thetaaT sigmaaT thetabT sigmabT thetarT sigmarT taur0 taur1 thrT sgmrt phirT ...
    from to

n=init(1);
h=init(2);
rf=init(3);
rs=init(4);
ca=init(5);
v=init(6);
na=init(7);
e=init(8);
rT=init(9);
hp=init(10);

% Na+ and K+ Equations and Currents
minf = 1/(1+exp((v-thetaM)/sigmaM));
ninf = 1/(1+exp((v-thetaN)/sigmaN));
tauN = tauNbar./cosh((v-thetaN)/(2*sigmaN));

alphaH = 0.128*exp(-(v+50)/18);
betaH = 4/(1+exp(-(v+27)/5));
hinf = alphaH/(alphaH+betaH);

iNa = gNa*(minf^3)*h*(v-vNa);
iK = gK*(n^4)*(v-vK);

% L-Type Ca++ Equations and Current
sinf = 1/(1+exp((v-thetaS)/sigmaS));
iCa = gCa*(sinf^2)*v*(Ca_ex/(1-exp((2*v)/RTF)));

% T-Type Ca++ Equations and Current
aTinf = 1/(1+exp((v-thetaaT)/sigmaaT));
bTinf = 1/(1+exp((rT-thetabT)/sigmabT))-1/(1+exp(-thetabT/sigmabT));
rTinf = 1/(1+exp((v-thetarT)/sigmarT));
taurT = taur0+taur1/(1+exp((v-thrT)/sgmrt));
iT=gT*(aTinf^3)*(bTinf^3)*v*(Ca_ex/(1-exp((2*v)/RTF)));


% SK Equations and Current
kinf = (ca^2)/(ca^2+(ks^2));
iSK = gSK*kinf*(v-vK);

% Na+ Persistant Current and Equations
mpinf = 1/(1+exp((v-thetaMP)/sigmaMP));
hpinf = 1/(1+exp((v-thetaHP)/sigmaHP));
tauhp = tauHPbar./cosh((v-thetaHP)/(2*sigmaHP));
iNap = gNap*mpinf*hp*(v-vNa);

% Na+ dependent K+ Current and Equations
winf = pmax/(1+(ec/na)^nh);
phiNa = (na^3)/(na^3+kp^3);
phiNaEqui = (neq^3)/(neq^3+kp^3);

iKNa = gKNa*winf*(v-vK);

% Hyperpolarization activated inward current Ih
rinff = 1./(1+exp((v-thetaRF)/sigmaRF));
rinfs = 1./(1+exp((v-thetaRS)/sigmaRS));
taurf = prf./(-7.4*(v+70)./(exp(-(v+70)/0.8)-1)+65*exp(-(v+56)/23));


iH = gH*(kr*rf+(1-kr)*rs)*(v-vH);


% A-type potassium current Ia
ainf = 1/(1+exp((v-thetaA)/sigmaA));
einf = 1/(1+exp((v-thetaE)/sigmaE));

iA = gA*ainf*e*(v-vK);

% Leak current
iL = gL*(v-vL);


if(time>=from&&time<=to)
    iap=iapp;
else
    iap=0;
end

output(1,1)=(ninf-n)/tauN;
output(2,1)=(hinf-h)/tauH;
output(3,1)=(rinff-rf)/taurf;
output(4,1)=(rinfs-rs)/taurs;
output(5,1)=-f*(eps*(iCa+iT)+ kca*(ca-0.1));
output(6,1)=(-iNa-iK-iCa-iT-iSK-iNap-iKNa-iH-iA-iL+iap)/C;
output(7,1)=-alphaNa*(iNa+iNap)-3*rpump*(phiNa-phiNaEqui);
output(8,1)=(einf-e)/tauE;
output(9,1)=phirT*(rTinf-rT)/taurT;
output(10,1)=(hpinf-hp)/tauhp;

