(Array indices must be positive integers or logical values.) any help plz

clear
%Propriétés des matériaux 
R1=630;  %Rayon de la roue en mm 
R2=inf ;  % largeur du rail en mm 
E1=210000; %Module d’élasticité du matériau 1 en MPa 
E2=210000; %Module d’élasticité du matériau 2 en MPa
v1=0.3; %Coefficient de Poisson du matériau 1
v2=0.3; %Coefficient de Poisson du matériau 2
u=0.5 %Coefficient de frottement 
p= 104125; % Charge appliquée en N 
L=135; %largeur entre roue/rail en mm  
q=p/L ; %charge appliqué entre roue rail en N/mm
%Mode de la sortie 
%Input choice=1 pour résultat adimensionné 
%Choice=2 pour résultat normale (avec dimension)
choice=2;
R=1/(1/R1+1/R2); %Rayon composé en mm
E=1/((1-v1^2)/E1+(1-v2^2)/E2); %Module d’élasticité composé en MPa
%Rayon de contact
a=(4*q*R/(pi*E))^(1/2); %Rayon de contact en mm
x=[-2*a:.01*3*a:2*a]; %Discrétisation de l’axe x
z=[0:.005*3*a:2*a]; % Discrétisation de l’axe z
p0=(q*E/(pi*R))^(1/2); %Pression maximal pour contact linéique en MPa 
%Réparation de la pression.
P(i)=p0*[1-(x./a).^2].^(1/2);
% Ci-dessous la boucle qui calcule les matrices de contrainte sous la surface :
s=[1/2*[((a^2)-(z.^2)-((x.^2))).^2+4*(a^2)*(x.^2).^1/2-((a^2)-(x.^2)-((z.^2)))]].^(1/2);
sx=(u*p0*(2*x./a)).*(1-(s./((a^2+s.^2).^(1/2))))+((x.*z.^2.*s.*a)./((a^2+s.^2).^(1/2)).*(s.^4+z.^2.*a^2)) ;
sy=(-2*v1*x./a).*(1-(s./((a^2+s.^2)))) ;
sz=(-u*p0)+((x.*z.^2.*s.*a)./ ((a^2+s.^2).^(1/2)).*(s.^4+a^2*z.^2)) ;
sxz=((u*p0)*(z./a)).*(2-(s./((a^2+s.^2).^(1/2))- ((a^2+s.^2).^(1/2)./s).*(x.^2.*s.^3*a.^2)./(a^2+s.^2).^(3/2).*(s.^4+z.^2*a^2))) ;
%OUTPUT
if choice==1
figure('name',' Contours contraintes sous surface xz, Modèle de Hertz');
subplot(231)
contour(xx/a, -zz/a, sx/p0)
C = contour(xx/a, -zz/a, sx/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('Contrainte Sx')
subplot(232)
contour(xx/a, -zz/a, sy/p0)
C = contour(xx/a, -zz/a, sy/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('Contrainte Sy')
subplot(233)
contour(xx/a, -zz/a, sz/p0)
C = contour(xx/a, -zz/a, sz/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('Contrainte Sz')
subplot(234)
contour(xx/a, -zz/a, sxz/p0)
C = contour(xx/a, -zz/a, sxz/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('txz')
subplot(235)
contour(xx/a, -zz/a, sxz/p0)
C = contour(xx/a, -zz/a, sxz/p0);
clabel(C)
xlabel('x/a')
ylabel('z/a')
title('tmax')
subplot(236)
plot(x/a,P./p0)
title('Pression normale')
xlabel('x/a')
ylabel('P/p0')
end
if choice==2
figure('name',' Contours contraintes sous surface xz: Modèle de Hertz');
subplot(231)
contour(xx, -zz, sx)
C = contour(xx, -zz, sx);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('Sxx [MPa]')
subplot(232)
contour(xx, -zz, sy)
C = contour(xx, -zz, sy);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('Syy [MPa]')
subplot(233)
contour(xx, -zz, sz)
C = contour(xx, -zz, sz);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('Szz [MPa]')
subplot(234)
contour(xx, -zz, txz)
C = contour(xx, -zz, sxz);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('txz [MPa]')
subplot(235)
contour(xx, -zz, sxz)
C = contour(xx, -zz, sxz);
clabel(C)
xlabel('x en mm')
ylabel('z en mm')
title('tmax [MPa]')
subplot(236)
plot(x,P)
title('P[MPa]')
xlabel('x en mm')
ylabel('P en MPa')
end

Answer 1

[the cyclist]

The specific reason you get that error is that sx is complex number.

sx is complex because s is complex.

s is complex because when you take the square root, elements 66 and 67 are negative, so you are taking the square root of a negative number.

I did not try to figure out any more than that.

If you don't know some of the tools for debugging MATLAB, you might want to read this documentation.

[Mokhles Elouaer]

thank you so much