(*M_P=1*)

(*Chaotic Potential m=10^-3*)

V1[ϕ_]=0.5*10^(-6) ϕ^2;                      
Plot[10^(6)V1[ϕ],{ϕ,-15,15},AxesLabel->{ϕ,"10^(6)""V[ϕ]"}, PlotStyle->{Blue,Thickness[0.005]},AxesStyle->{{Arrowheads[{0.,0.03}],{Thickness[0.005]}},{Arrowheads[{0.0,0.03}],{Thickness[0.005]}}},LabelStyle->Black,TicksStyle->Black]

(*Condiciones iniciales*)

a0chaotic=1; 
tlimchaotic=300000;

(*Solución para todas las condiciones iniciales posibles mediante NDSolve. (Flat Universe)*)

solutionchaotic=ParametricNDSolve[{ϕ''[t]+3(a'[t]/a[t]) ϕ'[t]+D[V1[ϕ[t]],ϕ[t]]==0,(a'[t]/a[t])==Sqrt[0.5ϕ'[t]^2+V1[ϕ[t]]],ϕ[0]==phi0chaotic,ϕ'[0]==phidot0chaotic,a[0]==a0chaotic},{ϕ,a},{t,0,tlimchaotic},{phi0chaotic,phidot0chaotic}];

(*Plot de la solución para el campo*)
Plot[{ϕ[15,0][t]/.solutionchaotic,0.1},{t,0,45000},PlotRange->All,PlotStyle->{{Green,Thickness[0.005]},{Black,Thickness[0.002]}},FrameLabel->{t,ϕ},Frame->True,FrameStyle->{{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black}}]

(*Plot logarítmico la solución para el factor de escala*)
LogPlot[{a[15,0][t]/.solutionchaotic,0.1},{t,0,45000},PlotRange->All,PlotStyle->{{Red,Thickness[0.006]},{Black,Thickness[0.002]}},FrameLabel->{t,a},Frame->True,FrameStyle->{{Thickness[0.008],Black},{Thickness[0.008],Black},{Thickness[0.008],Black},{Thickness[0.008],Black}}]

(*Cálculo del número de e-folds*)
asolutionchaotic[t_]=a[15,0][t]/.solutionchaotic;
Hchaotic[t_]=asolutionchaotic'[t]/asolutionchaotic[t];
NIntegrate[Hchaotic[t],{t,0,32100}]

(*--------------------------------------------------------------------------------------*)

(*Starobinski Potential m=10^(-6)*)

V2[ϕ_]=10^(-6)(-1+Exp[Sqrt[2/3]ϕ])^2;
Plot[10^(6)V2[ϕ],{ϕ,-10,1},AxesLabel->{ϕ,"10^(6)""V[ϕ]"}, PlotStyle->{Blue,Thickness[0.005]},AxesStyle->{{Arrowheads[{0.,0.03}],{Thickness[0.005]}},{Arrowheads[{0.0,0.03}],{Thickness[0.005]}}},LabelStyle->Black,TicksStyle->Black]

(*Condiciones Iniciales*)
a0starob=1; 
tlimstarob=2000000;

(*Solución para todas las condiciones iniciales posibles mediante ParametricNDSolve. (Flat Universe)*)

solutionstarob=ParametricNDSolve[{ϕ''[t]+3(a'[t]/a[t]) ϕ'[t]+D[V2[ϕ[t]],ϕ[t]]==0,(a'[t]/a[t])==Sqrt[0.5ϕ'[t]^2+V2[ϕ[t]]],ϕ[0]==phi0starob,ϕ'[0]==phidot0starob,a[0]==a0starob},{ϕ,a},{t,0,tlimstarob},{phi0starob,phidot0starob}];

(*Plot de la solución para el campo*)
Plot[{ϕ[-5,0][t]/.solutionstarob,0.00001},{t,0,150000},PlotRange->All,PlotStyle->{{Green,Thickness[0.006]},{Black,Thickness[0.002]}},FrameLabel->{t,ϕ},Frame->True,FrameStyle->{{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black}}]

(*Plot logarítmico la solución para el factor de escala*)
LogPlot[a[-5,0][t]/.solutionstarob,{t,0,150000},PlotRange->All,PlotStyle->{{Red,Thickness[0.006]},{Black,Thickness[0.002]}},FrameLabel->{t,a},Frame->True,FrameStyle->{{Thickness[0.008],Black},{Thickness[0.008],Black},{Thickness[0.008],Black},{Thickness[0.008],Black}}]

(*Cálculo del número de e-folds*)
asolution[t_]=a[-5,0][t]/.solutionstarob;
H[t_]=asolution'[t]/asolution[t];
NIntegrate[H[t],{t,0,134000}]