% State variables
v=x(1,1);  % voltage
n=x(2,1);  % open fraction
c=x(3,1);  % cytosol Calcium
cER = x(4,1); % ER Calcium

vk=-75; vl=-75; vca=25;
i=0; gk=2700; gl=150; gca=1000; cm=5300;
v1=-20.0; v2=24; v3=-16; v4=11.2; phi=0.035;

gkca=2000.0; 
Kca=5.0;

minf=.5*(1+tanh((v-v1)/v2));
ninf=.5*(1+tanh((v-v3)/v4));
tau=1/(cosh((v-v3)/(2*v4)));


alpha=4.50e-6; 
 
ikca =gkca/(1+(Kca/c))*(v-vk);
ica = gca*minf*(v-vca);

P_IP3R = 0.0008;
v_LSP = 0.6;
sigma = 0.032;
lambdaER = 2;
vlpm = 0.13;
f= 0.01;

% Vector field
xdot(1,1)=(+gl*(vl-v)+gk*n*(vk-v)- ica - ikca)/cm;
xdot(2,1)=phi*(ninf-n)/tau;
xdot(3,1)=f*(-alpha*ica - vlpm*c) + f/lambdaER*(P_IP3R*(100*cER-c)-v_LSP*c );    
xdot(4,1)=-1/100*f/(sigma*lambdaER)*(P_IP3R*(100*cER-c)-v_LSP*c );