(*Data generation for Fig3.12*)

Clear["Global`*"]

(*Load *)

P1=8*EE*II*(2*II/Ac)^(1/2)/l^3*u^3*(3/2-1/2*Tanh[u]^2-3/2*Tanh[u]/u)^(-1/2)

(* Deflection *)

delta1=P1*l^3/(4*EE*II)*(u-Tanh[u])/u^3

(*Axial Force*)

S1=(Ac*EE*P1^2/2*(3/2-1/2*Tanh[u]^2-3/2*Tanh[u]/u))^(1/3)

(* Deflection on linear theory *)

deltalin=P1*l^3/(12*EE*II)

(*Stress p.87 *)

fmax=1/3*EE*hl^2*u^2*(1+6^(1/2)*Tanh[u]*(3/2-1/2*Tanh[u]^2-3/2*Tanh[u]/u)^(-1/2))

flin=P1*l/(2*b*h*h/6)

(*specific data,unit [N], [mm] *)

EE=206*10^3; II=0.06937;Ac=3.226;l=67.2;hl=0.51/67.2;b=6.35;h=0.51;

P1

delta1

lstpd={};lstpdlin={};lstpS={};lstpfmax={};
Print["Load    Def.(lin.)    Def.(non-lin.)    Axial F."];
Do[
lstpd=Append[lstpd,{N[P1],delta1}];
lstpdlin=Append[lstpdlin,{N[P1],deltalin}];
lstpS=Append[lstpS,{N[P1],S1}];
lstpfmax=Append[lstpfmax,{N[P1],fmax}];
Print[N[P1]," , ",deltalin," , ",delta1," , ",S1],{u,0.001,2,0.1}];
ListLinePlot[{lstpd,lstpdlin},Frame->True,GridLines->Automatic,FrameLabel->{"Load","Deflection"},PlotStyle->{Thick,Dashed},LabelStyle->Directive[Red,Bold]]

ListLinePlot[lstpS,Frame->True,GridLines->Automatic,FrameLabel->{"Load","Tension"},LabelStyle->Directive[Red,Bold]]

lstpfmax={};lstpflin={};
Print["Load  Stress fmax(lin.)  Stress fmax(non-lin.) "];
Do[
lstpflin=Append[lstpflin,{N[P1],flin}];lstpfmax=Append[lstpfmax,{N[P1],fmax}];
Print[N[P1],"  ,  ",flin,"  ,  ",fmax ],{u,0.001,2,0.1}];
ListLinePlot[{lstpfmax,lstpflin},Frame->True,GridLines->Automatic,PlotStyle->{Thick,Dashed},FrameLabel->{"Load","Stress"},LabelStyle->Directive[Red,Bold]]