Why I am getting this message for Luo-Rudy Model?

Hi, I was trying to solve the Luo-Rudy model.These are the equations of the model. I used Euler method for the same.

This is the program I wrote for the same.

I got error like this.

%Luo-Rudy model using Euler method
clc; clear all; close all;
V_init=-84.5286;    %y(1)=V  y(2)=m  y(3)=h y(4)=j y(5)=d y(6)=f y(7)=X y(8)=Cai
m_init=0.0017;
h_init=0.9832;
J_init=0.995484;
d_init=0.000003;
f_init=1;
X_init=0.0057;
Cai_init=0.00002;
[V,m,h,J,d,f,X,Cai,t]=LRrun(0.1,500,V_init,m_init,h_init,J_init,d_init,f_init,X_init,Cai_init);
function [V,m,h,J,d,f,X,Cai,t]=LRrun(I,tspan,vi,mi,hi,Ji,di,fi,Xi,Cai_i)
%I=stim1(t(1),amp,n,a,c,J,T);
dt=0.001;
loop=ceil(tspan/dt);
EK=-77.5552;
Gbar_K=0.282;
Gbar_K1=0.6047;
EK1=-87.8789;
EKp=EK1;
%Initializing variable vectors
t=(1:loop)*dt;
V=zeros(loop,1);
m=zeros(loop,1);
h=zeros(loop,1);
J=zeros(loop,1);
d=zeros(loop,1);
f=zeros(loop,1);
X=zeros(loop,1);
Cai=zeros(loop,1);
%Set initial values for variables
V(1)=vi;
m(1)=mi;
h(1)=hi;
J(1)=Ji;
d(1)=di;
f(1)=fi;
X(1)=Xi;
Cai(1)=Cai_i;
%Euler method
for i=1:loop-1
    V(i+1)=V(i)+dt*(-(INa(V(i),m(i),h(i),J(i)) ...
        +Isi(V(i),d(i),f(i),Cai(i))+Gbar_K*X(i)*Xi(V)*(V(i)-EK)+Gbar_K1*(alphaK1(V(i))/(alphaK1(V(i)) ...
        +betaK1(V(i))))*(V(i)-EK1)+0.0183*Kp(V(i))*(V(i)-EKp)+0.03921*(V(i)+59.87))+I);
    m((i+1))=m(i)+dt*(alphaM(V(i))*(1-m(i))-betaM(V(i))*m(i)); 
    h((i+1))=h(i)+dt*(alphaH(V(i))*(1-h(i))-betaH(V(i))*h(i));
    J((i+1))=J(i)+dt*(alphaJ(V(i))*(1-J(i))-betaJ(V(i))*J(i));
    d((i+1))=d(i)+dt*(alphad(V(i))*(1-d(i))-betad(V(i))*d(i));
    f((i+1))=f(i)+dt*(alphaf(V(i))*(1-f(i))-betaf(V(i))*f(i));
    X((i+1))=X(i)+dt*(alphaX(V(i))*(1-X(i))-betaX(V(i))*X(i));
    Cai((i+1))=10^-4*(V(i)-(7.7-13.0287*log(Cai(i))))+0.07*(10^-4-Cai(i));
end
end
function ina=INa(V,m,h,J)
Gbar_Na=23;
E_Na=50;
ina=Gbar_Na*m^3*h*J*(V-E_Na);
end
function isi=Isi(V,d,f,Cai)
tol=10^-10;
isi=0.09*d*f*(V-(7.7-13.02878*log(max(Cai,tol))));
end
function xi=Xi(V)
 if V<=-100
    xi=1;  
 else
    xi=2.837*(exp(0.04*(V+77))-1)/(V+77)*exp(0.04*(V+35));  
 end
end
function aK1=alphaK1(V)
EK1=-87.8789;
aK1=1.02/(1+exp(0.2835*(V-EK1-59.215)));
end
function bK1=betaK1(V)
EK1=-87.8789;
bK1=(0.49124*exp(0.080328*(V-EK1+5.476))+(exp(0.061758*(V-EK1-594.31))))/(1+exp(-0.5143*(V-EK1+4.753)));
end
function kp=Kp(V)
kp=1/(1+exp(1+exp(7.488-V)/5.98));
end
function aM=alphaM(V)
aM=0.32*(V+47.13)/(1-exp(-0.1*(V+47.13)));
end
function bM=betaM(V)
bM=0.08*exp(-V/11);
end
function aH=alphaH(V)
if V>=-40
    aH=0;
else
    aH=0.135*exp((80+V)/(-6.8)); 
end
end
function bH=betaH(V)
if V>=-40
    bH=1/(0.13*(1+exp((V+10.66)/-11.1)));
else
    bH=3.56*exp(0.079*V)+3.1*10^5*exp(0.35*V);
end
end
function aJ=alphaJ(V)
if V>=-40
    aJ=0;
else
    aJ=(-1.2714*10^5*exp(0.2444*V)-3.474*10^-5*exp(-0.04391*V))*(V+37.78)/(1+exp(0.311*(V+79.23)));
end
end
function bJ=betaJ(V)
if V>=-40
    bJ=0.3*exp(-2.535*10^-7*V)/(1+exp(-0.1*(V+32)));
else
    bJ=0.1212*exp(-0.01052*V)/(1+exp(-0.1378*(V+40.14)));
end
end
function ad=alphad(V)
ad=0.095*exp(-0.01*(V-5))/(1+exp(-0.072*(V-5)));
end
function bd=betad(V)
bd=0.07*exp(-0.017*(V+44))/(1+exp(0.05*(V+44)));
end
function af=alphaf(V)
af=0.012*exp(-0.008*(V+28))/(1+exp(0.15*(V+28)));
end
function bf=betaf(V)
bf=0.0065*exp(-0.02*(V+30))/(1+exp(-0.2*(V+30)));
end
function aX=alphaX(V)
aX=0.0005*exp(0.083*(V+50))/(1+exp(0.057*(V+50)));
end
function bX=betaX(V)
bX=0.0013*exp(-0.06*(V+20))/(1+exp(-0.04*(V+20)));
end

Array indices must be positive integers or logical values.

Error in LReuler>LRrun (line 45)
    V(i+1)=V(i)+dt*(-(INa(V(i),m(i),h(i),J(i))+Isi(V(i),d(i),f(i),Cai(i))+Gbar_K*X(i)*Xi(V)*(V(i)-EK)+Gbar_K1*(alphaK1(V(i))/(alphaK1(V(i))+betaK1(V(i))))*(V(i)-EK1)+0.0183*Kp(V(i))*(V(i)-EKp)+0.03921*(V(i)+59.87)));

Error in LReuler (line 11)
[V,m,h,J,d,f,X,Cai,t]=LRrun(0.1,500,V_init,m_init,h_init,J_init,d_init,f_init,X_init,Cai_init);
 
>> 

Answer 1

[Alan Stevens]

You have a clashing definition for Xi. Once as a constant and once as a function. I suggest you change the function name to, say Xii and change the call in the following line

+Isi(V(i),d(i),f(i),Cai(i))+Gbar_K*X(i)*Xi(V)*(V(i)-EK)+Gbar_K1*(alphaK1(V(i))/(alphaK1(V(i))

to

+Isi(V(i),d(i),f(i),Cai(i))+Gbar_K*X(i)*Xii(V)*(V(i)-EK)+Gbar_K1*(alphaK1(V(i))/(alphaK1(V(i))