mistake ( Z must be at least a 2x2 matrix.)

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 p= 104125; % Charge appliquée en N

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

u=0.5 ; %coefficient de frottement choisit

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

a=(4*q*R/(pi*E)).^(1/2) ; %Rayon de contact en mm

x=[-3*a:.01*3*a:3*a]; %Discrétisation de l’axe x

z=[0:.005*3*a:3*a]; % Discrétisation de l’axe

p0=(q*E/(pi*R)).^(1/2); %Pression maximal pour contact linéique en MPa.

Pr=p0*(1-(x/a).^2).^0.5; %Réparation de la pression

sxz = ((u*p0)/a).*(-a+(3/2).*z.*(cot(a./z))-a.*(z.^(2))./(2*(a^2.+z.^2)));

%OUTPUT

figure('name',' Contours contraintes sous surface xz, Modèle de Hertz');

plot(sxz)

contour(x/a,-z/a,sxz/p0)

C = contour(x/a,-z/a,sxz/p0);

clabel(C)

xlabel('x/a')

ylabel('z/a')

title('Contrainte sxz')

grid on ;

figure('name','Contours contraintes sous surface xz: Modèle de Hertz');

plot(sxz)

contour(x,-z,sxz)

C = contour(x,-z,sxz);

clabel(C)

xlabel('x en mm')

ylabel('z en mm')

title('sxz [MPa]')

grid on ;

[Mokhles Elouaer]

any help please