Problème d'indice - equation de la chaleur - scilab - différences finies

[invitec2c8051e]

Bonjour, je fais actuellement un petit programme de résolution de la chaleur sous scilab par la méthode des différences finies.
J'ai pas mal debuggé mon code jusqu'à maintenant, mais il me donne en erreur que mes tailles de matrice sont mauvaises, j'en suppose que je suis bloqué pour une question d'indice puisque normalement mes matrices sont à la bonne tailles. je ne sais pas si je me suis trompé dans la structure de l'algorithme (Gérer les boucles) mais je pense que ça peut être ça.

Je n'ai pas vraiment de connaissance en programmation, ça doit être un des premiers que je fais et encore moins en méthode des différences finies, donc j'ai peut être manqué quelque chose de simple. Je commence à m'arracher les cheveux pour cet indice (for m=1 à la ligne 36) et surement les suivants...

Si vous pouviez m'aider ça serait super cool
Merci beaucoup.

function [eqchal]=simulachaleur(Cp1, Rho1, K1, h, Cp2, Rho2, K2, Tcons, Tini, V, Nbtube, etube, eair, temps, dt, dx)
funcprot(0)
//création matrice vide
T=[];
//nombre d'itération
boucle=ceil(temps/dt);
//nombre noeud total
noeud=ceil((etube+eair)/dx);
//nombre de noeud avant la limite
limite=ceil(etube/dx);
disp(boucle)
//Condition initial
for m=1:noeud
T(1,m)=Tini;
end

//Tracer les graphes/preparation
show_window(); clf()
f=gcf();
f.color_map=jetcolormap(256);
drawlater();
colorbar(Tini,Tcons);

// recopie la matrice vide avec condition initial
T_new=T; 

//Facteur de stabilité
alpha=K1/(Rho1*Cp1)
lambda=dt/(dx*dx)
if lambda<(1/(2*alpha)) then disp("Méthode instable, rentrer un autre dx et dt")
    break
end
//Méthode des différences finis
for n=1:boucle
    //Condition à gauche
    for m=1
        T_new(n+1,m)=(T(n,2)/dt)+((Tcons*(h/K1))/(h/K1+1/dt))
    end
    //Conduction dans la matière jusqu'à H2
	for m=2:limite
        T_new(n+1,m)=T(n,m)+K1*dt*(T(n,m+1)+T(n,m-1)-2*T(n,m))/(dt*dt*Rho1*Cp1)
    end
    //Conservation du flux
	for m=limite+1
		T_new(n+1,m)=((-K1/dx)*T(n,m-1)+(K2/dx)*T(n,m+1))
    end
    //Conduction dans l'H2
	for m=limite+1:noeud-1
	T_new(n+1,m)=T(n,m)+K2*dt*(T(n,m+1)+T(n,m-1)-2*T(n,m))/(dt*dt*Rho2*Cp2)
    end
    //Flux nul
	for m=noeud
	T_new(n+1,m)=T(n,m-1)
	end

//création des axes
x=1:noeud;
y=1:boucle;
z=T_new();
disp(z)
//Trace le graphique
Sgrayplot(x,y,z,strf="042",zminmax=[Tini,Tcons])
xtitle("simulation de la chaleur")
drawnow();
end
endfunction

-----