Programmation scilab des systémes déquations

Bonjour,
Voici que mon encadreur ma dit

ton code ne marche toujours pas...Il faut éviter d'utiliser des indices nuls de tableaux (en scilab, un tableau commence toujours avec un indice=1). Mais ce n'est pas le seul bug, notamment tu dois aussi bien initialiser et bien utiliser le tableau des x.

//Données numériques
N=11; // nombre de discrétisations en N pts donc (N-1) intervalles
nt=10;//nombre de temps effectués
tf=20;// temps final
t0=0;// temps initial
x0=0; //position initiale
xF=10;// position finale
deltat=(tf-t0)/nt;// pas de temps
deltax=(xF-x0)/(N-1);// pas d'espace


//Données biologiques
mu=1.;
tau=0.1;
E=15.;
R=0.2;
D=0.01;
s=1.;
eps=0.01;
r=0.5;

//Constructions de la matrice M
M=zeros(N,N);
for i=1:N-1
M(i,i)=2;
M(i,i+1)=-1;
M(i+1,i)=-1;
end

M(N,N)=2;

//Fatorisation de la matrice M en Lu
[L,U]=lu(M);

//Construction de la fonction Heaviside
function [H]=Heaviside(x)
if x>1
H(x-1)=1
elseif x<1
H(x-1)=0
else
H(x-1)=0.5
end
endfunction

//Construction de u,n et rho, à l'instant KT=0

rho_i=zeros(N-1,1);
if rho_i>0 & rho_i<1 then
rho(x,0)= rho_i+(1-rho_i)*H(x-1);
end
//********** Boucle en temps*********////
// initialisation
u=zeros(N,N);
rho=zeros(N,N);
n=zeros(N,N);
if i==0
n(0,k)=n(1,k);
rho(0,k)=rho(1,k);
u(0,k)=0;
end
if i==N
u(N,k)=0;
rho(N,k)=rho(N-1,k);
n(N,k)=n(N-1,k);
end

//k; itération en temps courante


for k=1:nt+1
for i=2:N-3

//print(%io(2),u,rho,n);


//print(%io(2),i,k);
u(i,k+1)= deltax^2*deltat*1/mu*s*rho(i,k)*u(i,k)-tau*deltax^2/mu*(n(i+1,k)*rho(i+1,k)/R^2+rho(i,k)^2-n(i,k)*rho(i,k)/R^2+rho(i,k)^2)+...
(1-E*deltat/mu)*(u(i,k)-2*u(i,k)+u(i-1,k));
u(N)=0;

rho(i,k+1)=rho(i,k)+eps.*n(i,k)*(1-rho(i,k))-1/deltax*(rho(i+1,k)*(u(i+1,k+1)-u(i+1,k))-rho(i,k)*(u(i,k+1)-u(i,k)));
rho(N-1)=rho(N);

n(i,k+1)=n(i,k)-1/deltax*(n(i+1,k)*(u(i+1,k+1)-u(i+1,k))-n(i,k)*(u(i,k+1)-u(i,k)))-deltat/deltax^2*(n(i+1,k)*(rho(i+2,k)-...
rho(i+1,k)-n(i,k)*(rho(i+1,k)-rho(i,k))))+deltat/deltax^2*D*(n(i+2,k)-2*n(i+1,k)+n(i,k))+deltat*r*n(i,k)*(1-n(i,k));
n(N-1)=n(N);

end

// Construction de F le second membre
F=zeros(N,N);
for i=2:N-3

// print(%io(2),F);

F(i,k)=deltax^2*deltat*1/mu*s*rho(i,k)*u(i,k)-tau*deltax^2/mu*(n(i+1,k)*rho(i+1,k)/R^2+rho(i,k)^2-n(i,k)*rho(i,k)/R^2+rho(i,k)^2)+...
(1-E*deltat/mu)*(u(i,k)-2*u(i,k)+u(i-1,k));
end
//mise a jour de la matrice M
M=zeros(N,N);
// print(%io(2),M);
for i=1:N-1
M(i,i)=2;
M(i,i+1)=-1;
M(i+1,i)=-1;
end

M(N,N)=2;
// mise a jour de la factorisation LU

[L,U]=lu(M);

// résolution de l'équation Mu=F
// ulocal la solution du systéme d'équation Mu=F
ulocal =zeros(N,1);
//print(%io(2),ulocal)
ulocal=M\F;
//mise a jour de rho et n

for i=2:N-3
rho(i,k+1)=rho(i,k)+eps.*n(i,k)*(1-rho(i,k))-1/deltax*(rho(i+1,k)*(u(i+1,k+1)-u(i+1,k))-rho(i,k)*(u(i,k+1)-u(i,k)));
n(i,k+1)=n(i,k)-1/deltax*(n(i+1,k)*(u(i+1,k+1)-u(i+1,k))-n(i,k)*(u(i,k+1)-u(i,k)))-deltat/deltax^2*(n(i+1,k)*(rho(i+2,k)-rho(i+1,k)-n(i,k)*(rho(i+1,k)-rho(i,k))))+deltat/deltax^2*D*(n(i+2,k)-2*n(i+1,k)+n(i,k))+deltat*r*n(i,k)*(1-n(i,k));

end
end