(*Modelo Aurrekoetxea, M_P=1*)
ClearAll["Global`*"]

(*Definimos el potencial y las ecuaciones del instanton a resolver*)
j=1;
k=1.5;
l=10^(-3);

V[ϕ_]=l (ϕ^2-6^2+ϕ^5)^2+k/Cosh[ϕ+j]^20;


EQ1[x_]=ϕ''[x]+3 (r'[x]/r[x]) ϕ'[x]-V'[ϕ[x]];
EQ2[x_]=(r''[x]/r[x])+(1/3)*(V[ϕ[x]]+ϕ'[x]^2);

ϕ[x_]=ϕ0+ϕ02*x^2;
r[x_]=x+r03*x^3;


ϕ2=ϕ02/.Solve[FullSimplify[Series[EQ1[x],{x,0,1}]]==0,ϕ02];
r3=r03/.Solve[FullSimplify[Series[EQ2[x],{x,0,1}]]==0,r03];

ϕini[x_]=ϕ0+ϕ2*x^2;
rini[x_]=x+r3*x^3;

(*Condiciones de contorno*)
xi=1/10^3;
xf=10;

Clear[ϕ,r]

BC={ϕ[xi]==ϕini[xi],ϕ'[xi]==ϕini'[xi],r[xi]==rini[xi],r'[xi]==rini'[xi]};

(*Resolvemos para todos los valores de ϕ0*)Solution=ParametricNDSolve[{EQ1[x]==0,EQ2[x]==0,BC},{ϕ,r},{x,xi,xf},{ϕ0}];

(*Plot de la solución del instanton con el ϕ0 correcto*)
Plot[ϕ[-0.4279703625555217605391788993074442259967329][x]/.Solution,{x,xi,3.50375},PlotRange->All,PlotStyle->{Green,Thickness[0.005]},FrameLabel->{ξ,ϕ},Frame->True,FrameStyle->{{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black}}]

(*Plot de la solución del factor de escala del exterior con el ϕ0 correcto*)
Plot[r[-0.4279703625555217605391788993074442259967329][x]/.Solution,{x,xi,3.50375},PlotRange->All,PlotStyle->{Red,Thickness[0.005]},FrameLabel->{ξ,ρ},Frame->True,FrameStyle->{{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black}}]

(*Plot del potencial con los puntos del tunneling*)
Show[Plot[{V[ϕ]},{ϕ,-3,3},PlotRange->{0,6},AxesLabel->{ϕ,"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],ListPlot[{{-0.4279703625555217605391788993074442259967329,V[-0.4279703625555217605391788993074442259967329]},{-1.46026,V[-1.46026]}},PlotStyle->{Red,Red}]]

--------------------------------------------------------------------------------------------------------------------------

(*Obtenemos la evolución del campo y el factor de escala en el interior de la burbuja*)
(*Initial conditions*)
ti=10^(-7);

a0open=ti;(*Expansión de Taylor para t<<, a va como t*)
tlimopen=600;

(*Equation of motion*)solutionopen=ParametricNDSolve[{ϕ''[t]+3 (a'[t]/a[t]) ϕ'[t]+D[V[ϕ[t]],ϕ[t]]==0,(a'[t]/a[t])==Sqrt[0.5 ϕ'[t]^2+V[ϕ[t]]+a[t]^(-2)],ϕ[0]==ϕ0open,ϕ'[0]==ϕdot0open,a[0]==a0open},{ϕ,a},{t,ti,tlimopen},{ϕ0open,ϕdot0open}];

(*Evolución del campo en el interior de la burbuja*)
Plot[{ϕ[-0.42797036255552176053917889999,0][t]/.solutionopen},{t,ti,tlimopen},PlotRange->All,PlotStyle->{{Green,Thickness[0.005]}},FrameLabel->{t,ϕ},Frame->True,FrameStyle->{{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black}}]

(*Plot logarítmico de la solución del factor de escala del interior de la burbuja*)
LogPlot[{a[-0.42797036255552176053917889999,0][t]/.solutionopen,0.000001},{t,ti,tlimopen},PlotRange->All,PlotStyle->{{Red,Thickness[0.006]},{Black,Thickness[0.002]}},FrameLabel->{t,a},Frame->True,FrameStyle->{{Thickness[0.006],Black},{Thickness[0.006],Black},{Thickness[0.006],Black},{Thickness[0.006],Black}}]

(*Obtención del número de e-folds*)
asolutionopen[t_]=a[-0.1,0][t]/.solutionopen;
Hopen[t_]=asolutionopen'[t]/asolutionopen[t];
efoldsopen=NIntegrate[Hopen[t],{t,0,110}]