For[i=-5,i<2,i++,
r0=1;
m=2500;
μ0=0;
l=20*r0;
μ[z_]=μ0+m*z;
ν=0.49;
P0=100;
M0=0;
ρs=1;
ρp=10;
λ=2*ν*μ[z]/(1-2*ν);
Ip=π/4*r0^4;
ω=0;
γ1new=1;
γ1old=0.001;
γ4new=1;
γ4old=0.001;
rm=15*r0;

λμ=2*ν/(1-2*ν);

kr=10^(i);

Ep=kr/(Ip/(μ[l]*2*(1+ν)*l^4));


While[Abs[(γ1old-γ1new)/γ1old]>0.001&&Abs[(γ4old-γ4new)/γ4old]>0.001,

γ1=γ1new;
γ2=(λμ+1)/(λμ+2);
γ3=(λμ+3)/(λμ+2);

γ4=γ4new;
γ5=(λμ+1);
γ6=(λμ+3);

s=First[NDSolve[{r*ϕ1''[r]+ϕ1'[r]-(γ1/r0)^2*r*ϕ1[r]-γ2*ϕ2'[r]-γ3*(ϕ1[r]-ϕ2[r])/r==0,r*ϕ2''[r]+ϕ2'[r]-
(γ4/r0)^2*r*ϕ2[r]+γ5*ϕ1'[r]+γ6*(ϕ1[r]-ϕ2[r])/r==0,ϕ1[r0]==1,ϕ2[r0]==1,ϕ1[rm]==0,ϕ2[rm]==0},
{ϕ1,ϕ2},{r,r0,rm},MaxSteps->∞]];

m1=π*ρs*Integrate[(((ϕ1[r]/.s)^2+(ϕ2[r]/.s)^2)*r),{r,r0,rm},Assumptions->r0>0]+π*r0^2*ρp//N;

c=Pi*Integrate[((ϕ1[r]/.s)^2+(ϕ2[r]/.s)^2)*r,{r,r0,rm},Assumptions->r0>0]//N;

k=π*(λμ+2)*Integrate[(((∂_r) (ϕ1[r]/.s))^2*r),{r,r0,rm},
Assumptions->r0>0]+π*(λμ+3)*Integrate[r*(((( ϕ1[r]- ϕ2[r])/r))^2)/.s,{r,r0,rm},
Assumptions->r0>0]+
π*Integrate[(((∂_r)( ϕ2[r]/.s))^2*r),{r,r0,rm},Assumptions->r0>0]+2*π*(λμ)*
Integrate[((ϕ1[r]- ϕ2[r])*(∂_r) ϕ1[r])/.s,{r,r0,rm},Assumptions->r0>0]+2*π*
Integrate[((ϕ1[r]- ϕ2[r])*(∂_r) ϕ2[r])/.s,{r,r0,rm},Assumptions->r0>0]//N;

m2=π*r0^2*ρs+π*ρs*Integrate[(((ϕ1[r]/.s)^2+(ϕ2[r]/.s)^2)*r),{r,r0,rm},Assumptions->r0>0]//N;

cs=Pi*Integrate[(((ϕ1[r]/.s)^2+(ϕ2[r]/.s)^2)*r),{r,r0,rm},Assumptions->r0>0]+π*r0^2//N;
h=cs*m;

u=First[y/.NDSolve[{Ep*Ip*y''''[z]-c*μ[z]*y''[z]+(k*μ[z]-m1*ω^2)*y[z]==0,
Ep*Ip*y'''[0]-c*μ[0]*y'[0]-P0==0,y''[l]==0,Ep*Ip*y''[0]+M0==0,Ep*Ip*y'''[l]-c*μ[l]*y'[l]+
(1/(Sqrt[cs]))(Sqrt[k]) *μ[l] (-cs-((h-(((Sqrt[cs]) m2 (ω^2))/(Sqrt[k]))) HypergeometricU[1+((h-(((Sqrt[cs]) m2 (ω^2))/(Sqrt[k])))/(2 cs m)),2,((2 (Sqrt[k]) *μ[l])/((Sqrt[cs]) m))])/
(m HypergeometricU[((h-(((Sqrt[cs]) m2 (ω^2))/(Sqrt[k])))/(2 cs m)),1,((2 (Sqrt[k])*μ[l])/((Sqrt[cs]) m))])) y[l]==0},y,{z,0,l},MaxSteps->∞]];

;
u2[z_]=(((E^(((Sqrt[k]) (l-z))/(Sqrt[cs]))) HypergeometricU[((h-(((Sqrt[cs]) m2 (ω^2))/(Sqrt[k])))/(2 cs m)),1,((2 (Sqrt[k]) μ[z])/((Sqrt[cs]) m))] u[l])/(HypergeometricU[((h-(((Sqrt[cs]) m2 (ω^2))/(Sqrt[k])))/(2 cs m)),1,((2 (Sqrt[k]) μ[l])/((Sqrt[cs]) m))]));

n1=Integrate[2*μ[z]*((∂_z) u[z])^2,{z,0,l}]+Integrate[2*μ[z]*((∂_z) u2[z])^2,{z,l,∞}]//N;
n2=Integrate[2*ρs*ω^2*(u[z])^2,{z,0,l}]+Integrate[2*ρs*ω^2*(u2[z])^2,{z,l,∞}]//N;

n3=Integrate[(λ+2*μ[z])*(u[z])^2,{z,0,l}]+Integrate[(λ+2*μ[z])*(u2[z])^2,{z,l,∞}]//N;
n4=n3/((λμ+2));
γ1new=(Sqrt[(n1+n2)/n3*r0^2]);
γ4new=(Sqrt[(n1+n2)/n4*r0^2]);
γ1old=γ1;
γ4old=γ4

];
Umax[i]=u[0];
Ih[i]=u[0]/(P0/(μ[l]*2*(1+ν)*l));
kri[i]=kr;
Mmax[i]=First[FindMaximum[{Ep*Ip*u''[z],0≤z≤l}, {z,l/2}]]
]