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

(*Definimos el potencial y las ecuaciones diferenciales y resolvemos los desarrollos de Taylor.*)
V[ϕ_,b_]=-1/2 ϕ^2-0.6/3 ϕ^3+1/4 ϕ^4+b;


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

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

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

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

(*Plot de la solución para el campo*)
Plot[{ϕ[0.7,0.34016506231][x]/.Solution},{x,xi,6.47},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 para el factor de escala en el exterior*)
Plot[{r[0.7,0.34016506231][x]/.Solution},{x,xi,6.47},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[ϕ,0.7]},{ϕ,-2,2},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.3401604307206375`,V[0.3401604307206375`,0.7]},{-0.29935276041568304`,V[-0.29935276041568304`,0.7]}},PlotStyle->{Red,Red}]]

(*Factor de escala y evolución del campo en el interior*)
(*Initial conditions*)
ti=10^(-7);
a0inside=ti;
tliminside=20;

(*Equation of motion*)

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

Plot[{ϕin[0.3401604307206375`,0][t]/.solutioninside},{t,ti,tliminside},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[{ain[0.3401604307206375`,0][t]/.solutioninside},{t,ti,tliminside},PlotRange->All,PlotStyle->{{Red,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}}]