(*******************  filename: IterateS.m (S := R/Rmin)  ********************)

<<DifferentialEquations`InterpolatingFunctionAnatomy`

targettime = targettimesec = 2 10^20

targettimeyr = targettimesec/(60*60*24*365)

Print[StringForm[" "]]

Print[StringForm["targettime in seconds = ``;  targettime in years = ``",
                  N[N[targettimesec,1000]]  ,  N[N[targettimeyr,1000]]   ]]


NDSolve`Iterate[state, targettime]

sol = NDSolve`ProcessSolutions[state]

         SSol = S /. sol

      SdotSol = Sdot /. sol

    Ssol[t_] := S[t] /. sol

    Rsol[t_] := N[Rmin Ssol[t],1000]

 Sdotsol[t_] := Sdot[t] /. sol

   T1sol[t_] := N[N[USimplify[c^2 A/3 Second^2],1000]] Ssol[t]

   T2sol[t_] := (N[USimplify[2 c^2 B Second^2],1000]/
                             N[(3 Rmin^6),1000])/Ssol[t]^5


Sddotsol[t_] := T1sol[t] + T2sol[t]

 Hsqrsol[t_] := USimplify[c^2 (N[A/3,1000] - N[k/Rsol[t]^2,1000]
                                    - N[B/(3 Rsol[t]^6),1000])]

    Hsol[t_] := USimplify[Sqrt[Hsqrsol[t]]]

    Qsol[t_] := (Sddotsol[t]/Ssol[t])/(Sdotsol[t]/Ssol[t])^2


{lowtime, hightime} = InterpolatingFunctionDomain[First[S /. sol]][[1,1]]


Print[StringForm[" "]]

Print[StringForm["lowtime = ``,  hightime = ``",
                  N[lowtime]  ,  N[hightime]    ]]

Print[StringForm[" "]]

Print[StringForm["Ssol[targettime] = ``,  Rsol[targettime] = ``",
                  N[Ssol[targettime]]  ,  N[Rsol[targettime]]    ]]

Print[StringForm[" "]]

Print[StringForm["Hsol[targettime] = `` = ``",
                  N[Hsol[targettime]],
                  N[Convert[Hsol[targettime],
                            (Kilo Meter)/(Mega Parsec Second)]] ]]


Print[StringForm[" "]]

Print[StringForm["Memory in use = ``",
                  N[MemoryInUse[]/2^20,7] megabytes]]

(*************************  End of file IterateS.m  **************************)
