Scilab résolution d'équation différentielle

Bonjour,

voici le programme que j'ai élaboré pour résoudre 2 équations différentielles par la méthode RK4 : (pour ne pas encombrer j'ai pas mis tous les paramètres). le probleme réside essentiellement dans les 2 dernières lignes.merci de m'aider.

function f=fun (y,t)
f(1)=-k*2*Concb^2*(1-y(1))*(M-y(1)); // la derivée par rapport à la fraction molaire
f(2)=((-ra*V*DHr)-Ua*(y(2)-Tp)-Qevac)/(CpA*Noncb+CpB*Nnh3+Nh2o*CpW); // la derivée par rapport à la temperature
endfunction

for i=1:1:n // la grande boucle : variation du temps
for j=1:1:nb // K1 pour les deux équations
ya(j)=y(i,j);
end
K1=fun(ya,t(i));

for j=1:1:nb //
k2(j) = Dt * K1(j);
yb(j)= y(i)+k2(j);
t(i)=t(i)+Dt;
end
K2=fun(yb,t(i));

for j=1:1:nb
k3(j) = Dt * K2(j);
yc(j)= y(i)+k3(j);
t(i)=t(i)+Dt;
end
K3=fun(yc,t(i));

for j=1:1:nb
k4(j)=h*K3(j);
yd(j)=y(i)+k4(j);
t(i)=t(i)+h;
end
K4=fun(yd,t(i));

p=(K1+2*K2+2*K3+K4)/6;

y(i+1)=y(i)+p(1,i)*h; //le probleme est dans cette ligne

t(i+1)=t(i)+h;

end