Clear[b,p0]
$Assumptions={Element[x|z|xx|zz,Reals], b>0,p0>0};
bxz=b^2-x^2+z^2;
xz=4*x^2*z^2;
koren = Sqrt[bxz^2+xz];
m2=1/2*(koren+bxz);
n2=1/2*(koren-bxz);
m=Sqrt[m2];
n=Sign[x]Sqrt[n2];
mpn=m2+n2;
sxx=-p0/b*(m*(1+((z^2+n2)/mpn))+2*z)// Simplify;
syy=-p0/b*m*(1-((z^2+n2)/mpn))// Simplify;
sxy=p0/b*n*((m2-z^2)/mpn) // Simplify;
sigma={{sxx,sxy},{sxy,syy}};
psd=Sqrt[(((sxx-syy)/2))^2+sxy^2]/.Sign[x]^2->1;

(* izrazi za resitev v midplanu posebej iz knjige *)
y= z/b;
symxx=-p0((1+2y^2)/Sqrt[(1+y^2)]+2y);
symzz=-p0 /Sqrt[1+y^2];
tau1sym=-p0(y+y^2/Sqrt[1+y^2]);

(* podatki *)
h = 0.01;
w = 2h;

rl= {p0 -> 8.36 10^7,b -> 4.13 10^(-4)};

(* div sigma == 0 *)
dd1=Div[sigma[[1]],{x,z}];
dd2=Div[sigma[[2]],{x,z}];
plotdd1=dd1/. rl/.Derivative[1][Sign][x]->0;
plotdd2=dd2/. rl/.Derivative[1][Sign][x]->0;
plots={
Plot3D[plotdd1,{x,-w,w},{z,-h,h},ImageSize->Medium],
Plot3D[plotdd2,{x,-w,w},{z,-h,h},ImageSize->Medium]
};
Grid[{plots}]
plot=Plot3D[{plotdd1,plotdd2},{x,-w,w},{z,-h,h},ImageSize->Large, AxesLabel->{"x","y","divσ"}]
Export["images/divs_hertzian_analytical.png",plot];

(* Checks on midplane: general == special *)
(syy/.{x->0})  == symzz // FullSimplify (* mora biti true *)
zacyy=syy/.{z->0} // Simplify (* robni pogoj je eliptičen *)
Plot[zacyy/.rl,{x,-3b/.rl,3b/.rl}]
sxy/.{z->0} //Simplify(* mora biti 0 *)
sxy/.{x->0} (* mora biti 0 *)
psymxx=sxx/.{x->0}//Simplify; (* za sxx pokazimo numericno ujemanje (Simplify ni dovolj pameten) *)
Plot[{psymxx/.rl,symxx/.rl},{z,-h,0},PlotStyle->{{Dashed,Red},{Dotted,Green}}]
Plot[Evaluate[{psd/.{x->0}/.rl,tau1sym/.rl}],{z,-h,0},PlotStyle->{{Dashed,Red},{Dotted,Green}}] (* za psd tudi numericno *)

(* narisimo vse resitve v okolici pritiska normalizirane na p0 *)
plyy = Simplify[syy/p0/.{x->xx b,z->zz b}]/.rl;
plxx=Simplify[sxx/p0/.{x->xx b,z->zz b}]/.rl;
plxy=Simplify[sxy/p0/.{x->xx b,z->zz b}]/.rl;
ppsd=Simplify[psd/p0/.{x->xx b,z->zz b}]/.rl;
plots = {{
ContourPlot[plxx,{xx,-3,3},{zz,0,-3},AspectRatio->Automatic,ImageSize->600,PlotRange->All,FrameLabel->{"x/b","y/b"},PlotLabel->"(Subscript[s, xx])/(Subscript[p, 0])",Contours->20],
ContourPlot[plyy,{xx,-3,3},{zz,0,-3},AspectRatio->Automatic,ImageSize->600,PlotRange->All,FrameLabel->{"x/b","y/b"},PlotLabel->"(Subscript[s, yy])/(Subscript[p, 0])",Contours->20],
ContourPlot[plxy,{xx,-3,3},{zz,0,-3},AspectRatio->Automatic,ImageSize->600,PlotRange->All,FrameLabel->{"x/b","y/b"},PlotLabel->"(Subscript[s, xy])/(Subscript[p, 0])",Contours->20]},{
ContourPlot[ppsd,{xx,-3,3},{zz,0,-3},AspectRatio->Automatic,ImageSize->600,PlotRange->All,FrameLabel->{"x/b","y/b"},PlotLabel->"(Subscript[τ, 1])/(Subscript[p, 0]) = principal stress difference = max shear stress.",Contours->20]
}};
Grid[plots]
plots = Flatten[plots];
Export["images/sxx_hertzian_analytical.png",plots[[1]]];
Export["images/syy_hertzian_analytical.png",plots[[2]]];
Export["images/sxy_hertzian_analytical.png",plots[[3]]];
Export["images/psd_hertzian_analytical.png",plots[[4]]];

pp0=p0/.rl;
plt = ParametricPlot[{{(plyy/.x->0)/pp0,z},{(plxx/.x->0)/pp0,z},{(ppsd/.x->0)/pp0,z},{Norm[sigma/.x->0/.rl]/pp0,z}},{z,-h,0},AspectRatio->1,Axes->False,GridLines->Automatic,Frame->True,PlotRange->All,PlotLegends->{"(Subscript[s, yy])","(Subscript[s, xx])","psd = (Subscript[τ, 1])","norm"},PlotLabel->"Stresses in the midplane x=0.", FrameLabel->{"stress / (Subscript[p, 0])","y"},ImageSize->Large]
Export["images/midplane_stress_hertzian_analytical.png",plt]