probleme avec MATLAB

[ABN84]

bonsoir,
en essayant d'executer le script ci dessous j'obtiens le message:

??? Subscript indices must either be real positive integers or logicals.

qqun saurait-il comment regler ça?
merci.

max=60;%nbr max d'iterations
l=6;%languer du mur
h=2;%hauteur du mur
e=0.5;%epaisseur du mur
n1=1/49;%n=1/(r-1)//r=50
n2=1/199;%n=1/(r-1)//r=200
n3=1/179;%n=1/(r-1)//r=200-20=180
lambdab=1.8;%lambde du beton
rho=2300;% rho du beton
c=800;
k=ones(50,200) * 273.13;%matrice qui sert à rammener les K en °C
sigma=5.67*10^(-8);
t0=298.13;
E1=1368;
%t1=323.13;
p0=t0^4*sigma;
s=l*e;
emis=0.9;

omega=ones(50,200) * 298.13;
d1=ones(50,1)*298.13;%base
d1(1,1)=298.13;
d1(50,1)=323.13;
d2=ones(1,200)*298.15;%est
d3=ones(50,1)*298.15;%haut
d4=ones(1,180)*298.15;%west
for m=1:1:max
        for i=2:1:49
            d1(i,1)=d1(i,1)+(-4*d1(i,1)+d1(i-1,1)+d1(i+1,1))/(2+n1*rho*c/lambdab);% conduction thermique sur la base
        end
end
d5=ones(1,20)*323.13;%west/radiateur T° cte au cours du temps
for a=0:1:180% a: angle que fait le soleil avec l'horizontale à partir de l'est
    if a<90% faces est  et superieure exposées au rayonnement
        t2=(((298.13^4)/2)+E1/(sigma*cos(a)))^0.25;% T° de la face est si on ne prends en compte que le rayonnement
        d2=ones(1,200)*t2;%est
        t3=(((298.13^4)/2)+E1/(sigma*cos(90-a)))^0.25;% T° de la face superieure si on ne prends en compte que le rayonnement
        d3=ones(50,1)*t3;%haut
        if a<45%rayonnement plus intance sur la face est que sur la face superieure=> la chaleur se deplace des points de la face est vers ceux de la face superieure => T° cte sur la face est correspondant au rayonnement + equation de la chaleur à une dimension sur la face superieure
           d3(1,1)=t2;
           for m=1:1:max
                for i=2:1:49
                   d3(i,1)=d3(i,1)+(-4*d3(i,1)+d3(i-1,1)+d3(i+1,1))/(2+n1*rho*c/lambdab);
                end
           end
        else%rayonnement plus intance sur la face sup que sur la face est => la chaleur se deplace des points de la face sup vers ceux de la face est => T° cte sur la face sup correspondant au rayonnement + equation de la chaleur à une dimension sur la face est           
           d2(1,1)=t3;
           for m=1:1:max
                for j=2:1:199
                   d2(1,j)=d2(1,j)+(-4*d2(1,j)+d2(1,j-1)+d2(1,j+1))/(2+n2*rho*c/lambdab);
                end
           end           
        end 
        d4(1,0)=d3(50,1);
        d4(1,180)=50;%face west non exposée au soleil
        for m=1:1:max
                for j=2:1:179
                   d4(1,j)=d4(1,j)+(-4*d4(1,j)+d4(1,j-1)+d4(1,j+1))/(2+n3*rho*c/lambdab);%eq de la chaleur à une dimension
                end
        end  
    else% faces sup et west exposées au rayonnement
        t4=(((298.13^4)/2)+E1/(sigma*cos(180-a)))^0.25;
        d4=ones(1,180)*t4;%west
        t3=(((298.13^4)/2)+E1/(sigma*cos(a-90)))^0.25;
        d3=ones(50,1)*t3;%haut
        if a<135
           d4(0,1)=t3;
           d4(180,1)=50;
           for m=1:1:max
                for j=2:1:179
                   d4(1,j)=d4(1,j)+(-4*d4(1,j)+d4(1,j-1)+d4(1,j+1))/(2+n3*rho*c/lambdab);
                end
           end    
        else
           d3(50,1)=t4;
           for m=1:1:max
                for i=2:1:49
                   d3(i,1)=d3(i,1)+(-4*d3(i,1)+d3(i-1,1)+d3(i+1,1))/(2+n1*rho*c/lambdab);
                end
           end 
        end
        d2(1,1)=d3(0,1);%face west non exposée au soleil
        d2(1,200)=d1(0,1);
           for m=1:1:max
                for j=2:1:199
                   d2(1,j)=d2(1,j)+(-4*d2(1,j)+d2(1,j-1)+d2(1,j+1))/(2+n2*rho*c/lambdab);%eq de la chaleur à une dimension
                end
           end 
        t2=(((298.13^4)/2)+E1/(sigma*cos(a)))^0.25;
        d2=ones(1,200)*t2;%est
    end
    d6=ones(1,200)*0;%d6 rassemble d4 et d5
    d6(1,1:20)=d5;
    d6(1,21:200)=d4;
    omega(50,:)=d6;
    omega(:,200)=d1;
    omega(1,:)=d2;
    omega(:,1)=d3;
    for j=2:1:199
        for i=2:1:49
            omega(i,j)=omega(i,j)+(-4*omega(i,j)+omega(i-1,j)+omega(i,j-1)+omega(i+1,j)+omega(i,j+1))/(4+n1*n2*rho*c/lambdab);%eq de la chaleur à deux dimensions
        end
    end    
    omegac=omega-k;
    %plot(omegac);
end

-----

[ABN84]

à mon sens, le probleme viendrait d'un conflit entre:
d4(0,1)=t3;
et
t3=(((298.13^4)/2)+E1/(sigma*cos(90-a)))^0.25;

en effet d4 requiere un reel, mais en affichant t3, je vois un complexe qui est dû au fait que cos(90-a)<0. et c'est ce qui est bizarre car 0<a<90 donc cos(90-a) ne devrait pas etre <0