(*Gravitational Instanton 1, M_P=1*)
ClearAll["Global`*"]

(*Definimos nuestras ecuaciones diferenciales y hacemos los desarrollos de Taylor de las soluciones*)

V[ϕ_,c_]=0.5 ϕ^2-0.5 ϕ^3+(c/8) ϕ^4;


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

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

ϕ2=ϕ02/.Solve[FullSimplify[Series[EQ1[x,c],{x,0,1}]]
==0,ϕ02];
r3=r03/.Solve[FullSimplify[Series[EQ2[x,c],{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]};

Plot[{V[ϕ,0.5],V[ϕ,0.6],V[ϕ,0.7]},{ϕ,-2,7},AxesLabel->{ϕ,"V[ϕ]"},PlotStyle->{{RGBColor[0,1,1],Thickness[0.005]},{RGBColor[0,0.5,1],Thickness[0.005]},{RGBColor[0,0,1],Thickness[0.005]}},AxesStyle->{{Arrowheads[{0.,0.03}],{Thickness[0.005]}},{Arrowheads[{0.0,0.03}],{Thickness[0.005]}}},LabelStyle->Black,TicksStyle->Thick]

(*Resolvemos la ecuación diferencial para todos los ϕ0*)
Solution=ParametricNDSolve[{EQ1[x,c]==0,EQ2[x,c]==0,BC},{ϕ,r},{x,xi,xf},{c,ϕ0},Method->{"Shooting","StartingInitialConditions"->{ϕ[xi]==ϕini[xi],ϕ'[xi]==D[ϕini[xi],x]}}];

(*Plot de los instantones gravitacionales*)
Plot[{ϕ[0.5,3.95][x]/.Solution,ϕ[0.6,3.81][x]/.Solution,ϕ[0.7,3.4103][x]/.Solution},{x,xi,12.5},PlotRange->All,AxesLabel->{ξ,ϕ},PlotStyle->{{RGBColor[0,0.6,0],Thickness[0.005]},{RGBColor[0.1,1,0.1],Thickness[0.005]},{RGBColor[0.6,1,0.6],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 factor de escala del exterior*)
Plot[{r[0.5,3.95][x]/.Solution,r[0.6,3.81][x]/.Solution,r[0.7,3.4103][x]/.Solution},{x,xi,12.5},PlotRange->All,AxesLabel->{ξ,ϕ},PlotStyle->{{RGBColor[0.6,0,0],Thickness[0.005]},{RGBColor[1,0.2,0.2],Thickness[0.005]},{RGBColor[1,0.6,0.6],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 de tunneling*)
Show[Plot[{V[ϕ,0.5],V[ϕ,0.6],V[ϕ,0.7]},{ϕ,-2,7},AxesLabel->{ϕ,"V[ϕ]"},PlotStyle->{{RGBColor[0,1,1],Thickness[0.005]},{RGBColor[0,0.5,1],Thickness[0.005]},{RGBColor[0,0,1],Thickness[0.005]}},AxesStyle->{{Arrowheads[{0.,0.03}],{Thickness[0.005]}},{Arrowheads[{0.0,0.03}],{Thickness[0.005]}}},LabelStyle->Black,TicksStyle->Black],ListPlot[{{3.9514424095676812,V[3.9514424095676812,0.5]},{0.0023541788099001653`,V[0.0023541788099001653,0.5]},{0.02304919700539819`,V[0.02304919700539819`,0.6]},{0.04976275735324205`,V[0.04976275735324205`,0.7]},{3.81,V[3.81,0.6]},{3.4103,V[3.41030,0.7]}},PlotStyle->{Red,Red}]]

(*Factor de escala en el interior de la burbuja*)
(*Condiciones iniciales*)
ti=10^(-7);
a0inside=ti;
tliminside=3.375;

(*Equation of motion*)
solutioninside=ParametricNDSolve[{ϕin''[t]+3 (ain'[t]/ain[t]) ϕin'[t]+D[V[ϕin[t],q],ϕin[t]]==0,(ain''[t]/ain[t])==-1/3*( ϕin'[t]^2-V[ϕin[t],q]),ϕin[ti]==ϕ0in,ϕin'[ti]==ϕdot0in,ain[ti]==a0inside,ain'[ti]==1},{ϕin,ain},{t,ti,tliminside},{ϕ0in,ϕdot0in,q}];

Plot[{ain[3.9514424095676812,0,0.5][t]/.solutioninside,ain[3.81,0,0.6][t]/.solutioninside,ain[3.4103,0,0.7][t]/.solutioninside},{t,ti,tliminside},PlotRange->{0,1.3},PlotStyle->{{RGBColor[0.6,0,0],Thickness[0.005]},{RGBColor[1,0.2,0.2],Thickness[0.005]},{RGBColor[1,0.6,0.6],Thickness[0.005]}},FrameLabel->{t,a},Frame->True,FrameStyle->{{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black},{Thickness[0.005],Black}}]