(*            function  f[l] = (l/lambda)^q             *)
(* l is the wing crack length *)
(* all quantities are dimensionless; no units *)

(* setting the constants *)
dim=2000; (* number of points in the lists; integer *)
lambda=10.; (* value of parameter lambda; real number *)
n=10; (* subcritical crack growth exponent; integer *)
dl=N[10/dim]; (* increment of wing crack length *)
(* in this example, the final wing crack length is 10 *)

(* calculating the lists (or vectors) of wing crack length, lVec *)
(* stress intensity factor, kVec, crack growth velocity, vVec, and time, tVec *)
lVec=Table[dl (i-1), {i,1,10 dim+1}];
kVec= Table[Sqrt[1.+(lVec[[i]]/lambda)^q]-Sqrt[lVec[[i]]],{i,1,10 dim+1}];
ij=1;
While[kVec[[ij]]>0,
nDim=ij;
ij+=1];
temp=Take[lVec,nDim];
lVec=temp;
temp=Take[kVec,nDim];
kVec=temp;
vVec=Table[kVec[[i]]^n,{i,nDim}];
tVec=Table[0.,{i,nDim}];
Do[
dt=NIntegrate[(Sqrt[1.+(l/m)^q]-Sqrt[l])^-n,{l,lVec[[ii]],lVec[[ii+1]]}];
tVec[[ii+1]]=tVec[[ii]]+dt,{ii,nDim-1}];
(* the stress intensity factor is not allowed to reach zero or become negative *)
tVec[[1]]=tVec[[2]]/100.;
(* the zero initial time is replaced by a small finite value *)

(* calculating the characteristic parameters of the log-log curves i figure 8 *)
dlvVec=Differences[Log[vVec]];
dltVec=Differences[Log[tVec]];
dx=1000;
nn=1;
While[dx>10.^-7,
dx=dltVec[[nn]];
nn+=1];
temp=Take[dlvVec,nn-1];
dlvVec=temp;
temp=Take[dltVec,nn-1];
dltVec=temp;
failureTime=tVec[[nn-1]];
Print[tVec[[nn]],"  ",tVec[[nn-2]],"  failure time     ",failureTime];
flag=0; (* flag is set at zero *)
Do[
m1=dlvVec[[i]]/dltVec[[i]];(* power law exponent point i *)
m2=dlvVec[[i+1]]/dltVec[[i+1]];(* power law at point i+1 *)
If[m2>m1&&flag==0,(* inflexion point is reached at first increase of power law exponent *)
inflexionTime=tVec[[i]];
Print[m1,"  ",m2,"   inflexion time    ",inflexionTime];
flag=1]; (* flag set at 1 when inflexion point is passed *)
If[m1 m2<0.,(* minimum crack velocity when power law exponent changes sign *)
minimumCrackVelocityTime=(tVec[[i]]+tVec[[i+1]])/2.;
Print[exp1,"  ",exp2,"    minimum velocity time     ",minimum CrackVelocityTime]];
If[dltVec[[i]]==dltVec[[i+1]]&&flag==1,
flag=0;
Break[]],{i,3,nn-2}];
Print["   ratio     ",failureTime/minimumCrackVelocityTime];