Dynamic AdS/CMT - Superconductor - Spectral - General Ansatz

Print["hmm"]

GRSetup2.nb - modified

(* Latest Version *)

metric={gtt,grr,gθθ,gϕϕ,gzz}

Gd=DiagonalMatrix[metric]

Gd[[1,2]]=grt

Gd[[2,1]]=grt

BuildMetric:=(Gu=FullSimplify[Inverse[Gd]];G={{u,d},IdentityMatrix[Nind]};)

BuildMetric

ind={t,r,θ,ϕ,z}

Nind=Length[ind]

getind[i1_]:=(If[IntegerQ[i1],i1,Position[ind,i1][[1,1]]])

lstind[i1_]:=(lsttmp2={};Do[AppendTo[lsttmp2,getind[i1[[lsttmp1,1]]]],{lsttmp1,1,Length[i1]}];lsttmp2)

up[i1_,i2_]:=(lsttmp2={};Do[AppendTo[lsttmp2,If[i1[[lsttmp1,2]]=!=i2[[lsttmp1]],If[i1[[lsttmp1,2]]===d,d,u],0]],{lsttmp1,1,Length[i1]}];lsttmp2)

MT[t1_,dim_,indpos_]:=(MTtmp1=0;If[Length[indpos]=!=dim,Print["Error in index position"]];Do[MTtmp1=Table[MTtmp1,{MTtmp3,1,Nind}],{MTtmp2,1,dim}];t1={indpos,MTtmp1})

ST[t1_,tind_,val_]:=(STtmp1=t1;If[up[tind,t1[[1]]]=!=Table[0,{Length[tind]}],Print["Error in Index Assignment!!"];];STtmp1[[2]]=ReplacePart[t1[[2]],val,lstind[tind]];STtmp1)

Edited the expression below as 7.0 seems to complain about Sum[x] whereas 6.0 would just return x (ie. if there is no sum index!)

GT[t1_,tind_]:=(GTtmp2=up[tind,t1[[1]]];GTtmp3="";GTtmp4="";GTtmp5="";GTtmp6=t1;
		Do[
			Switch[GTtmp2[[GTtmp1]],
				0,
				(GTtmp3=GTtmp3<>","<>ToString[getind[tind[[GTtmp1,1]]]];),
				u,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gu[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}"),
				d,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gd[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}";)];
			,{GTtmp1,1,Length[t1[[1]]]}];If[GTtmp4===""&&GTtmp5==="",ToExpression["GTtmp6[[2"<>GTtmp3<>"]]"],ToExpression["Sum[GTtmp6[[2"<>GTtmp3<>"]]"<>GTtmp4<>GTtmp5<>"]"]])

Trace, Outer Product and Covariant Derivative! - GREinstein.nb

Syntax:   trT[tensor,{{1,4},{2,3}}]   -  ie. contract 1 and 4 plus 2 and 3 positions

trT[t1_,lst_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[tracetmp];

If[dimt1===2,

tracetmp=Sum[GT[tmpt1,{{s1,d},{s1,u}}],{s1,1,Nind}];

,

MT[tracetmp,dimt1-2Length[lst],Table[d,{ii,1,dimt1-2Length[lst]}]];

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

If[{}===Position[lst,indno],
tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring3=tmpstring3<>"{i"<>ToString[indno]<>",1,Nind},";
];

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",d},";

,{indno,1,dimt1}];

tmpstring4="";

Do[

tmpstring4=tmpstring4<>"{s"<>ToString[indno]<>",1,Nind},";

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,1]]]->"s"<>ToString[indno]];

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,2]]]<>",d}"->"s"<>ToString[indno]<>",u}"];

,{indno,1,Length[lst]}];

tmpeval="Do[tracetmp=ST[tracetmp,{"<>StringDrop[tmpstring1,-1]<>"},Sum[GT[tmpt1,{"<>StringDrop[tmpstring2,-1]<>"}],"<>StringDrop[tmpstring4,-1]<>"]],"<>StringDrop[tmpstring3,-1]<>"]";

ToExpression[tmpeval];

];

tracetmp
)

Syntax:   outerT[tensor1,tensor2]

outerT[t1_,t2_]:=(

dimt1=Length[Dimensions[t1[[2]]]];dimt2=Length[Dimensions[t2[[2]]]];
tmpt1=t1;tmpt2=t2;
Clear[outertmp];

MT[outertmp,dimt1+dimt2,Table[d,{ii,1,dimt1+dimt2}]];

tmpstring1="";tmpstring2="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

tmpstring3="";tmpstring4="";

Do[

tmpstring3=tmpstring3<>"{j"<>ToString[indno]<>",d},";
tmpstring4=tmpstring4<>"{j"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt2}];

tmpeval="Do[outertmp=ST[outertmp,{"<>tmpstring1<>StringDrop[tmpstring3,-1]<>"},GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[tmpt2,{"<>StringDrop[tmpstring3,-1]<>"}]],"<>tmpstring2<>StringDrop[tmpstring4,-1]<>"]";

ToExpression[tmpeval];

outertmp

)

scalarDT[t1_]:=(

tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,1,{d}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{{i1,d}},D[tmpt1,ind[[i1]]]],{i1,1,Nind}]";

ToExpression[tmpeval];

covDTtmp

)

covDT[t1_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,dimt1+1,Table[d,{ii,1,dimt1+1}]];

(* Build expression *)

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

Do[
tmpstring4="-Sum[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[Γ,{{REPL,d},{s1,d},{s2,u}}],{s2,1,Nind}]";

tmpstring4=StringReplace[tmpstring4,"i"<>ToString[indno]-> "s2"];

tmpstring4=StringReplace[tmpstring4,"REPL"->"i"<>ToString[indno]];

tmpstring3=tmpstring3<>tmpstring4;

,{indno,1,dimt1}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{"<>tmpstring1<>"{s1,d}},D[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}],ind[[s1]]]"<>tmpstring3<>"],"<>tmpstring2<>"{s1,1,Nind}]";


ToExpression[tmpeval];

covDTtmp

);

.

Buildmetric

Gd=(1/(z^2))(({{-T[v,z], -gvz[v,z], 0, 0}, {-gvz[v,z], gzz[v,z], 0, 0}, {0, 0, S[v,z]^2, 0}, {0, 0, 0, S[v,z]^2}}))

ind={v,z,x1,x2}

Nind=Length[ind]

BuildMetric;

Gd

Gu

detg=Simplify[Det[Gd]]

GRrun2.nb

Clear[Γ]

MT[Γ,3,{d,d,u}];

Clear[R4]

MT[R4,4,{d,d,d,u}];

Clear[R2]

MT[R2,2,{d,d}];

Clear[G2]

MT[G2,2,{d,d}];

Do[Γ=ST[Γ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}},Sum[(1/2)GT[G,{{ind[[i3]],u},{ind[[ss1]],u}}](D[GT[G,{{ind[[i1]],d},{ind[[ss1]],d}}],ind[[i2]]]+D[GT[G,{{ind[[ss1]],d},{ind[[i2]],d}}],ind[[i1]]]-D[GT[G,{{ind[[i1]],d},{ind[[i2]],d}}],ind[[ss1]]]),{ss1,1,Nind}]],{i1,1,Nind},{i2,1,Nind},{i3,1,Nind}]

(* FullSimplify[Γ] *)

Γ=FullSimplify[Γ];

Do[
	R4=ST[R4,
			{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}},
	
			D[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i2]]]-D[GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i1]]]+
			Sum[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i2]],d},{ind[[ss1]],d},{ind[[i4]],u}}]-GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i1]],d},{ind[[ss1]],d},{ind[[i4]],u}}],{ss1,1,Nind}]
		
		]
		,
	{i1,1,Nind},{i2,1,Nind},{i3,1,Nind},{i4,1,Nind}]

(* FullSimplify[R4] *)

Do[R2=ST[R2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			Sum[GT[R4,{{ind[[i1]],d},{ind[[ss1]],u},{ind[[i2]],d},{ind[[ss1]],d}}],{ss1,1,Nind}]
		]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[R2] *)

Ricci=Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}];

Do[G2=ST[G2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			GT[R2,{{ind[[i1]],d},{ind[[i2]],d}}]-(1/2)GT[G,{{ind[[i1]],d},{ind[[i2]],d}}]Ricci]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[G2] *)

(* FullSimplify[Ricci] *)

Einstein equations

R2=Simplify[R2];

Ricci=Simplify[Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}]];

Clear[A];

MT[A,1,{d}];

A=ST[A,{{v,d}},av[v,z]];

A=ST[A,{{z,d}},az[v,z]];

Clear[F];

MT[F,2,{d,d}];

Do[
F=ST[F,{{i1,d},{i2,d}},D[GT[A,{{i1,d}}],ind[[i2]]]-D[GT[A,{{i2,d}}],ind[[i1]]]];
,{i1,1,Nind},{i2,1,Nind}];

trFF=Simplify[trT[outerT[F,F],{{2,4}}]];

Fsqr=Simplify[trT[trFF,{{1,2}}]];

ϕ[v_,z_]:=p[v,z]+I α[v,z]

Clear[Dp];

MT[Dp,1,{d}];

Do[
Dp=ST[Dp,{{i1,d}},D[ϕ[v,z],ind[[i1]]]-I GT[A,{{i1,d}}]ϕ[v,z]];
,{i1,1,Nind}];

ConjDp={{d},Conjugate[Dp[[2]]]};

DpDp=Simplify[ComplexExpand[trT[outerT[ConjDp,Dp],{1,2}]]];

Vterm=Simplify[-2ComplexExpand[Conjugate[ϕ[v,z]]ϕ[v,z]]];

DpDp=(DpDp/.{Re[xx_]->xx,Im[xx_]->0});

sqrtdetg=Simplify[PowerExpand[Sqrt[Det[-Gd]]]];

L=sqrtdetg(Ricci+6-(1/4)Fsqr-DpDp-Vterm)

tmp1=Simplify[Coefficient[L,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]];

tmp2=Simplify[Coefficient[L,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]];

tmp3=Simplify[Coefficient[L,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]];

tmp4=Simplify[Coefficient[L,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]];

GHBT=Simplify[D[tmp1.{Derivative[1,0][gvz][v,z],0,(1/2)Derivative[0,1][gvz][v,z]},v]+D[tmp1.{0,Derivative[0,1][gvz][v,z],(1/2)Derivative[1,0][gvz][v,z]},z]];

GHBT+=Simplify[D[tmp2.{Derivative[1,0][gzz][v,z],0,(1/2)Derivative[0,1][gzz][v,z]},v]+D[tmp2.{0,Derivative[0,1][gzz][v,z],(1/2)Derivative[1,0][gzz][v,z]},z]];

GHBT+=Simplify[D[tmp3.{Derivative[1,0][T][v,z],0,(1/2)Derivative[0,1][T][v,z]},v]+D[tmp3.{0,Derivative[0,1][T][v,z],(1/2)Derivative[1,0][T][v,z]},z]];

GHBT+=Simplify[D[tmp4.{Derivative[1,0][S][v,z],0,(1/2)Derivative[0,1][S][v,z]},v]+D[tmp4.{0,Derivative[0,1][S][v,z],(1/2)Derivative[1,0][S][v,z]},z]];

L2=L-GHBT;

tmp1=Simplify[Coefficient[L2,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]]

tmp2=Simplify[Coefficient[L2,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]]

tmp3=Simplify[Coefficient[L2,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]]

tmp4=Simplify[Coefficient[L2,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]]

EOM[fld_]:=(
tmp=D[L2,fld];
Do[
If[(tmp2=D[fld,ind[[ii]]])=!=0,
tmp-=D[D[L2,tmp2],ind[[ii]]];
]
,{ii,1,Nind}];
tmp
)

eqp=EOM[p[v,z]];

eqα=EOM[α[v,z]];

eqav=EOM[av[v,z]];

eqaz=EOM[az[v,z]];

eqgvz=EOM[gvz[v,z]];

eqgzz=EOM[gzz[v,z]];

eqT=EOM[T[v,z]];

eqS=EOM[S[v,z]];

Clear[gvz,gzz,az]

gvz[v_,z_]:=1

gzz[v_,z_]:=0

az[v_,z_]:=0

eqp2=Simplify[eqp]

eqα2=Simplify[eqα]

eqav2=Simplify[eqav]

eqaz2=Simplify[eqaz]

eqgvz2=Simplify[eqgvz]

eqgzz2=Simplify[eqgzz]

eqT2=Simplify[eqT]

eqS2=Simplify[eqS]

p[t_,r_]:=0

α[t_,r_]:=0

av[t_,r_]:=0

T[t_,r_]:=1

S[t_,r_]:=1

Simplify[eqp2]

Simplify[eqα2]

Simplify[eqav2]

Simplify[eqaz2]

Simplify[eqgvz2]

Simplify[eqgzz2]

Simplify[eqT2]

Simplify[eqS2]

Clear[α,p,av,T,S]

av[v_,z_]:=0

α[v_,z_]:=0

eqp3=Simplify[eqp]

eqgvz3=Simplify[eqgvz]

eqgzz3=Simplify[eqgzz]

eqT3=Simplify[eqT]

eqS3=Simplify[eqS]

Momsol=Simplify[Solve[{eqT3==0,eqgzz3==0},{Derivative[2,0][S][v,z], Derivative[0,2][S][v,z]}][[1]]]

Evsol=Simplify[Solve[{eqgvz3==0,eqS3==0,eqp3==0},{Derivative[1,1][p][v,z], Derivative[1,1][S][v,z], Derivative[0,2][T][v,z]}][[1]]]

Print["\n\nTO DO: check these agree with previous!\n\n"]

Clear[α,p,av,T,S]

Evsol=Simplify[Solve[{eqgvz2==0,eqS2==0,eqp2==0,eqα2==0,eqaz==0},{Derivative[1,1][p][v,z],Derivative[1,1][α][v,z], Derivative[1,1][av][v,z], Derivative[1,1][S][v,z], Derivative[0,2][T][v,z]}][[1]]]

Momsol=Simplify[Solve[{eqT2==0,eqgzz==0,eqav==0},{Derivative[2,0][S][v,z], Derivative[0,2][S][v,z], Derivative[0,2][av][v,z]}][[1]]]

Simplify[D[eqT2,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqav,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqgzz,z]/.D[Evsol,v]/.Evsol/.Momsol]

av[v_,z_]:=a[v,z]

α[v_,z_]:=h[v,z]

Const1=Derivative[0,2][S][v,z]-(Derivative[0,2][S][v,z]/.Momsol)

Const2=Derivative[0,2][a][v,z]-(Derivative[0,2][a][v,z]/.Momsol)

Const3=Derivative[2,0][S][v,z]-(Derivative[2,0][S][v,z]/.Momsol)

Eq1=Derivative[1,1][p][v,z]-(Derivative[1,1][p][v,z]/.Evsol)

Eq2=Derivative[1,1][h][v,z]-(Derivative[1,1][h][v,z]/.Evsol)

Eq3=Derivative[1,1][a][v,z]-(Derivative[1,1][a][v,z]/.Evsol)

Eq4=Derivative[1,1][S][v,z]-(Derivative[1,1][S][v,z]/.Evsol)

Eq5=Derivative[0,2][T][v,z]-(Derivative[0,2][T][v,z]/.Evsol)

Clear[deriv];

deriv[X_]:={X[v,z]->X,D[X[v,z],v]->ToExpression[ToString[X]<>"v"],D[X[v,z],z]->ToExpression[ToString[X]<>"v"],D[D[X[v,z],v],z]->ToExpression[ToString[X]<>"vz"],D[X[v,z],{v,2}]->ToExpression[ToString[X]<>"vv"],D[X[v,z],{z,2}]->ToExpression[ToString[X]<>"zz"]};

deriv[j]

Clear[pow];

pow:={Power[X_,n_]->If[Hold[MemberQ[{p,a,h,T,S,pv,av,hv,Tv,Sv,az,hz,Tz,Sz,z},X]&&n>1],ToExpression[ToString[X]<>ToString[n]],Power[X,n]]}

pow2:={Power[z,n_]->If[Hold[n<0],ToExpression["invz"<>ToString[-n]],Power[z,n]]}

repl=Join[deriv[p],deriv[a],deriv[h],deriv[T],deriv[S]];

tmp=ToString[CForm[Release[Release[(Const1//.repl)/.pow]/.pow2]]];

C1str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Const2//.repl)/.pow]/.pow2]]];

C2str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Const3//.repl)/.pow]/.pow2]]];

C3str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Eq1//.repl)/.pow]/.pow2]]];

Eq1str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Eq2//.repl)/.pow]/.pow2]]];

Eq2str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Eq3//.repl)/.pow]/.pow2]]];

Eq3str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Eq4//.repl)/.pow]/.pow2]]];

Eq4str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Release[Release[(Eq5//.repl)/.pow]/.pow2]]];

Eq5str=StringReplace[tmp,"*"->".*"]

p[v_,z_]:=(z p1[v]+z^2 p2[v]+z^3 p3[z]+z^4 p4[z]+z^5 p5[z]);

h[v_,z_]:=(z h1[v]+z^2 h2[v]+z^3 h3[z]+z^4 h4[z]+z^5 h5[z]);

a[v_,z_]:=0;

T[v_,z_]:=1;

S[v_,z_]:=1;

Simplify[Series[Eq1,{z,0,3}]]

Simplify[Series[Eq2,{z,0,3}]]

p[v_,z_]:=(z p1[v]+z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z h1[v]+z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 a2[v]+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1;

S[v_,z_]:=1;

Simplify[Series[Eq1,{z,0,2}]]

Simplify[Series[Eq2,{z,0,2}]]

Simplify[Series[Eq3,{z,0,0}]]

Simplify[Series[Const2,{z,0,0}]]

Clear[p,h,a,T,S]

p[v_,z_]:=(z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2a2[v]+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1;

S[v_,z_]:=1;

tmp1=Simplify[Series[Eq1,{z,0,1}]]

tmp2=Simplify[Series[Eq2,{z,0,1}]]

tmp3=Simplify[Series[Eq3,{z,0,2}]]

tmp4=Simplify[Series[Const2,{z,0,2}]]

p[v_,z_]:=(z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 a2[v]+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1;

S[v_,z_]:=1;

tmp1=Simplify[Series[Eq1,{z,0,1}]]

tmp2=Simplify[Series[Eq2,{z,0,1}]]

tmp3=Simplify[Series[Eq3,{z,0,2}]]

tmp4=Simplify[Series[Const2,{z,0,2}]]

Solve[{tmp1,tmp2,tmp3,tmp4}==0,{p3[v],h3[v]}]

pos = 1;

p[v_,z_]:=(z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 a2[v]+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1+z^2T2[v]+z^3 T3[v]+z^4 T4[v]+z^5 T5[v];

S[v_,z_]:=1+z^2S2[v]+z^3 S3[v]+z^4 S4[v]+z^5 S5[v];

tmp1=Series[Const1,{z,0,pos-1}]

tmp2=Series[Const2,{z,0,pos-1}]

tmp3=Series[Const3,{z,0,pos-1}]

tmp4=Series[Eq1,{z,0,pos-1}]

tmp5=Series[Eq2,{z,0,pos-1}]

tmp6=Series[Eq3,{z,0,pos-1}]

tmp7=Series[Eq4,{z,0,pos-1}]

tmp8=Series[Eq5,{z,0,pos-1}]

p[v_,z_]:=(z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]+z^6p6[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]+z^6 h6[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^4 a4[v]+z^5 a5[v]+z^6 a6[v];

T[v_,z_]:=1+z^3 T3[v]+z^4 T4[v]+z^5 T5[v]+z^6 T6[v];

S[v_,z_]:=1+z^4 S4[v]+z^5 S5[v]+z^6 S6[v];

pos=2;

repltmp={a1'[v]->0};

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl=Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{h3[v],p3[v]}][[1]]

pos=3;

repltmp=Join[{a1'[v]->0},repl,D[repl,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl2=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T3'[v],a4[v],h4[v],p4[v],S4[v],T4[v]}][[1]]]

pos=4;

repltmp=Join[{a1'[v]->0},repl,D[repl,v],D[repl,{v,2}],repl2,D[repl2,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl3=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T5[v],a5[v],h5[v],p5[v],S5[v],T5[v]}][[1]]]

pos=5;

repltmp=Join[{a1'[v]->0},repl,D[repl,v],D[repl,{v,2}],D[repl,{v,3}],repl2,D[repl2,v],D[repl2,{v,2}],repl3,D[repl3,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl4=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T6[v],a6[v],h6[v],p6[v],S6[v],T6[v]}][[1]]]

Simplify[Series[p[v,z],{z,0,3}]/.repltmp]

Simplify[Series[h[v,z],{z,0,3}]/.repltmp]

Simplify[Series[a[v,z],{z,0,4}]/.repltmp]

Simplify[Series[T[v,z],{z,0,4}]/.repltmp]

Simplify[Series[S[v,z],{z,0,4}]/.repltmp]

Normal[Simplify[Series[p[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[h[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[a[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[T[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[S[v,z],{z,0,5}]/.repltmp]]

p[v_,z_]:=(z^2) p2[v]+(z^3) (a0[v] h2[v]+p2'[v])+(1/6) (z^4) (3 a1[v] h2[v]+4 (-((a0[v])^2) p2[v]+h2[v] a0'[v]+2 a0[v] h2'[v]+p2''[v]))+(1/3) (z^5) (-((a0[v])^3) h2[v]-p2[v] T3[v]+2 a1[v] h2'[v]+3 a0'[v] h2'[v]-3 ((a0[v])^2) p2'[v]+h2[v] a0''[v]+a0[v] (-2 a1[v] p2[v]-3 p2[v] a0'[v]+3 h2''[v])+Derivative[3][p2][v]);

h[v_,z_]:=(z^2) h2[v]+(z^3) (-a0[v] p2[v]+h2'[v])+(1/6) (z^4) (-4 ((a0[v])^2) h2[v]-3 a1[v] p2[v]-4 p2[v] a0'[v]-8 a0[v] p2'[v]+4 h2''[v])+(1/3) (z^5) (((a0[v])^3) p2[v]-h2[v] T3[v]-3 ((a0[v])^2) h2'[v]-2 a1[v] p2'[v]-3 a0'[v] p2'[v]-p2[v] a0''[v]-a0[v] (2 a1[v] h2[v]+3 (h2[v] a0'[v]+p2''[v]))+Derivative[3][h2][v]);

a[v_,z_]:=a0[v]+z a1[v]+(1/6) (z^4) (a0[v] (((h2[v])^2)+((p2[v])^2))-p2[v] h2'[v]+h2[v] p2'[v])+(1/30) (z^5) (5 a1[v] (((h2[v])^2)+((p2[v])^2))+4 (((h2[v])^2) a0'[v]+p2[v] (p2[v] a0'[v]+2 a0[v] p2'[v]-h2''[v])+h2[v] (2 a0[v] h2'[v]+p2''[v])));

T[v_,z_]:=1+(1/12) (z^4) (3 ((a1[v])^2)-4 (((h2[v])^2)+((p2[v])^2)))+(z^3) T3[v]+(1/3) (z^5) (-h2[v] h2'[v]-p2[v] p2'[v]);

S[v_,z_]:=1+(1/6) (z^4) (-((h2[v])^2)-((p2[v])^2))-(3/10) (z^5) (h2[v] h2'[v]+p2[v] p2'[v]);

pos=4;

tmp1=Simplify[Series[Const1,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5,{z,0,pos-1}]]

Clear[p,h,a,T,S]

pos = 1;

p[v_,z_]:=(z p1[v]+z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 a2[v]+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1+z^2T2[v]+z^3 T3[v]+z^4 T4[v]+z^5 T5[v];

S[v_,z_]:=1+z^2S2[v]+z^3 S3[v]+z^4 S4[v]+z^5 S5[v];

tmp1=Series[Const1,{z,0,pos-1}]

tmp2=Series[Const2,{z,0,pos-1}]

tmp3=Series[Const3,{z,0,pos-1}]

tmp4=Series[Eq1,{z,0,pos-1}]

tmp5=Series[Eq2,{z,0,pos-1}]

tmp6=Series[Eq3,{z,0,pos-1}]

tmp7=Series[Eq4,{z,0,pos-1}]

tmp8=Series[Eq5,{z,0,pos-1}]

repl=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{a1'[v],S2[v],T2[v],a2[v]}][[1]]]

pos = 1;

p[v_,z_]:=(z p1[v]+z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 (-h2[v] p1[v])+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1+z^2(-(1/2) ((p1[v])^2))+z^3 T3[v]+z^4 T4[v]+z^5 T5[v];

S[v_,z_]:=1+z^2(-(1/4) ((p1[v])^2))+z^3 S3[v]+z^4 S4[v]+z^5 S5[v];

pos=2;

repltmp={a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])};

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl=Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{h3[v],p3[v],S3[v],a3[v]}][[1]]

pos=3;

repltmp=Join[{a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])},repl,D[repl,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl2=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T3'[v],a4[v],h4[v],p4[v],S4[v],T4[v]}][[1]]]

pos=4;

repltmp=Join[{a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])},repl,D[repl,v],D[repl,{v,2}],repl2,D[repl2,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl3=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T5[v],a5[v],h5[v],p5[v],S5[v],T5[v]}][[1]]]

repltmp=Join[{a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])},repl,D[repl,v],D[repl,{v,2}],D[repl,{v,3}],repl2,D[repl2,v],D[repl2,{v,2}],repl3,D[repl3,v]];

Simplify[Series[p[v,z],{z,0,3}]/.repltmp]

Simplify[Series[h[v,z],{z,0,3}]/.repltmp]

Simplify[Series[a[v,z],{z,0,4}]/.repltmp]

Simplify[Series[T[v,z],{z,0,4}]/.repltmp]

Simplify[Series[S[v,z],{z,0,4}]/.repltmp]

Normal[Simplify[Series[p[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[h[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[a[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[T[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[S[v,z],{z,0,5}]/.repltmp]]

p[v_,z_]:=z p1[v]+(z^2) p2[v]+(z^3) (a0[v] h2[v]+(((p1[v])^3)/2)+p2'[v])+(1/6) (z^4) (a1[v] (3 h2[v]-2 a0[v] p1[v])-4 ((a0[v])^2) p2[v]+5 ((p1[v])^2) p2[v]-p1[v] T3[v]+4 h2[v] a0'[v]+8 a0[v] h2'[v]+4 ((p1[v])^2) p1'[v]+4 p2''[v])+(1/12) (z^5) (-4 ((a0[v])^3) h2[v]-3 ((a1[v])^2) p1[v]-8 ((h2[v])^2) p1[v]+4 ((p1[v])^5)+6 p1[v] ((p2[v])^2)-4 p2[v] T3[v]+12 a0'[v] h2'[v]+a1[v] (-2 p1[v] a0'[v]+8 h2'[v])+4 p1[v] p2[v] p1'[v]-T3[v] p1'[v]+8 p1[v] ((p1'[v])^2)+4 ((a0[v])^2) (2 ((p1[v])^3)-3 p2'[v])+13 ((p1[v])^2) p2'[v]+4 h2[v] a0''[v]+a0[v] (17 h2[v] ((p1[v])^2)-4 (3 p2[v] a0'[v]+a1[v] (2 p2[v]+p1'[v])-3 h2''[v]))+3 ((p1[v])^2) p1''[v]+4 Derivative[3][p2][v]);

h[v_,z_]:=(z^2) h2[v]+(z^3) (-(1/2) a1[v] p1[v]-a0[v] p2[v]+h2'[v])+(1/6) (z^4) (-4 ((a0[v])^2) h2[v]+9 h2[v] ((p1[v])^2)-3 a1[v] p2[v]-4 p2[v] a0'[v]-2 a1[v] p1'[v]+a0[v] (3 ((p1[v])^3)-8 p2'[v])+4 h2''[v])+(1/12) (z^5) (4 ((a0[v])^3) p2[v]+14 h2[v] p1[v] p2[v]-4 h2[v] T3[v]+3 ((p1[v])^3) a0'[v]+2 ((a0[v])^2) (a1[v] p1[v]-6 h2'[v])+16 ((p1[v])^2) h2'[v]+20 h2[v] p1[v] p1'[v]-12 a0'[v] p2'[v]-4 p2[v] a0''[v]-2 a1[v] (3 ((p1[v])^3)+4 p2'[v]+p1''[v])+a0[v] (-8 a1[v] h2[v]+p1[v] T3[v]+((p1[v])^2) (-4 p2[v]+9 p1'[v])-12 (h2[v] a0'[v]+p2''[v]))+4 Derivative[3][h2][v]);

a[v_,z_]:=a0[v]+z a1[v]-(z^2) h2[v] p1[v]+(1/6) (z^3) p1[v] (3 a1[v] p1[v]+4 a0[v] p2[v]-4 h2'[v])+(1/12) (z^4) (4 ((a0[v])^2) h2[v] p1[v]+a0[v] (2 ((h2[v])^2)-3 ((p1[v])^4)+2 ((p2[v])^2)+8 p1[v] p2'[v])+2 (2 p1[v] p2[v] a0'[v]-p2[v] h2'[v]+a1[v] p1[v] (3 p2[v]+p1'[v])+h2[v] (-5 ((p1[v])^3)+p2'[v])-2 p1[v] h2''[v]))+(1/30) (z^5) (-4 ((a0[v])^3) p1[v] p2[v]-24 h2[v] ((p1[v])^2) p2[v]+3 h2[v] p1[v] T3[v]+4 ((h2[v])^2) a0'[v]-3 ((p1[v])^4) a0'[v]+4 ((p2[v])^2) a0'[v]+12 ((a0[v])^2) p1[v] h2'[v]-19 ((p1[v])^3) h2'[v]-16 h2[v] ((p1[v])^2) p1'[v]+12 p1[v] a0'[v] p2'[v]+4 p1[v] p2[v] a0''[v]-4 p2[v] h2''[v]+a1[v] (5 ((h2[v])^2)+9 a0[v] h2[v] p1[v]-2 ((a0[v])^2) ((p1[v])^2)+10 ((p1[v])^4)+5 ((p2[v])^2)+2 p2[v] p1'[v]+11 p1[v] p2'[v]+2 p1[v] p1''[v])+4 h2[v] p2''[v]+a0[v] (-((p1[v])^2) T3[v]+((p1[v])^3) (4 p2[v]-9 p1'[v])+8 (h2[v] h2'[v]+p2[v] p2'[v])+12 p1[v] (h2[v] a0'[v]+p2''[v]))-4 p1[v] Derivative[3][h2][v]);

T[v_,z_]:=1-(1/2) (z^2) ((p1[v])^2)+(z^3) T3[v]+(1/12) (z^4) (3 ((a1[v])^2)-2 (2 ((h2[v])^2)+6 a0[v] h2[v] p1[v]+((p1[v])^4)+2 p2[v] (p2[v]-2 p1'[v])+2 p1[v] p2'[v]))+(1/6) (z^5) (a1[v] p1[v] (-3 h2[v]+2 a0[v] p1[v])+4 ((a0[v])^2) p1[v] p2[v]-2 ((p1[v])^3) p2[v]+((p1[v])^2) T3[v]-h2[v] p1[v] a0'[v]-2 h2[v] h2'[v]+a0[v] (-5 p1[v] h2'[v]+3 h2[v] p1'[v])-2 p2[v] p2'[v]+3 p1'[v] p2'[v]-p1[v] p2''[v]);

S[v_,z_]:=1-(1/4) (z^2) ((p1[v])^2)-(1/3) (z^3) p1[v] p2[v]+(1/96) (z^4) (-16 ((h2[v])^2)-24 a0[v] h2[v] p1[v]-11 ((p1[v])^4)-16 ((p2[v])^2)-24 p1[v] p2'[v])+(1/60) (z^5) (a1[v] p1[v] (3 h2[v]+4 a0[v] p1[v])+8 ((a0[v])^2) p1[v] p2[v]-17 ((p1[v])^3) p2[v]+2 ((p1[v])^2) T3[v]-8 h2[v] p1[v] a0'[v]-18 h2[v] h2'[v]-16 a0[v] p1[v] h2'[v]-8 ((p1[v])^3) p1'[v]-18 p2[v] p2'[v]-8 p1[v] p2''[v]);

pos=4;

tmp1=Simplify[Series[Const1,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5,{z,0,pos-1}]]

crepl={a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v]),T3'[v]->-p1[v] (a0[v] h2[v]+((a0[v])^2) p1[v]+p2'[v]-p1''[v])};

Simplify[tmp1]

Simplify[tmp2]

Simplify[tmp3/.crepl]

Simplify[tmp4/.crepl]

Simplify[tmp5/.crepl]

Simplify[tmp6/.crepl]

Simplify[tmp7]

Simplify[tmp8]

Dynamic AdS/CMT - Superconductor - Spectral - Numerial Ansatz

Print["hmm"]

GRSetup2.nb - modified

(* Latest Version *)

metric={gtt,grr,gθθ,gϕϕ,gzz}

Gd=DiagonalMatrix[metric]

Gd[[1,2]]=grt

Gd[[2,1]]=grt

BuildMetric:=(Gu=FullSimplify[Inverse[Gd]];G={{u,d},IdentityMatrix[Nind]};)

BuildMetric

ind={t,r,θ,ϕ,z}

Nind=Length[ind]

getind[i1_]:=(If[IntegerQ[i1],i1,Position[ind,i1][[1,1]]])

lstind[i1_]:=(lsttmp2={};Do[AppendTo[lsttmp2,getind[i1[[lsttmp1,1]]]],{lsttmp1,1,Length[i1]}];lsttmp2)

up[i1_,i2_]:=(lsttmp2={};Do[AppendTo[lsttmp2,If[i1[[lsttmp1,2]]=!=i2[[lsttmp1]],If[i1[[lsttmp1,2]]===d,d,u],0]],{lsttmp1,1,Length[i1]}];lsttmp2)

MT[t1_,dim_,indpos_]:=(MTtmp1=0;If[Length[indpos]=!=dim,Print["Error in index position"]];Do[MTtmp1=Table[MTtmp1,{MTtmp3,1,Nind}],{MTtmp2,1,dim}];t1={indpos,MTtmp1})

ST[t1_,tind_,val_]:=(STtmp1=t1;If[up[tind,t1[[1]]]=!=Table[0,{Length[tind]}],Print["Error in Index Assignment!!"];];STtmp1[[2]]=ReplacePart[t1[[2]],val,lstind[tind]];STtmp1)

Edited the expression below as 7.0 seems to complain about Sum[x] whereas 6.0 would just return x (ie. if there is no sum index!)

GT[t1_,tind_]:=(GTtmp2=up[tind,t1[[1]]];GTtmp3="";GTtmp4="";GTtmp5="";GTtmp6=t1;
		Do[
			Switch[GTtmp2[[GTtmp1]],
				0,
				(GTtmp3=GTtmp3<>","<>ToString[getind[tind[[GTtmp1,1]]]];),
				u,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gu[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}"),
				d,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gd[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}";)];
			,{GTtmp1,1,Length[t1[[1]]]}];If[GTtmp4===""&&GTtmp5==="",ToExpression["GTtmp6[[2"<>GTtmp3<>"]]"],ToExpression["Sum[GTtmp6[[2"<>GTtmp3<>"]]"<>GTtmp4<>GTtmp5<>"]"]])

Trace, Outer Product and Covariant Derivative! - GREinstein.nb

Syntax:   trT[tensor,{{1,4},{2,3}}]   -  ie. contract 1 and 4 plus 2 and 3 positions

trT[t1_,lst_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[tracetmp];

If[dimt1===2,

tracetmp=Sum[GT[tmpt1,{{s1,d},{s1,u}}],{s1,1,Nind}];

,

MT[tracetmp,dimt1-2Length[lst],Table[d,{ii,1,dimt1-2Length[lst]}]];

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

If[{}===Position[lst,indno],
tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring3=tmpstring3<>"{i"<>ToString[indno]<>",1,Nind},";
];

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",d},";

,{indno,1,dimt1}];

tmpstring4="";

Do[

tmpstring4=tmpstring4<>"{s"<>ToString[indno]<>",1,Nind},";

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,1]]]->"s"<>ToString[indno]];

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,2]]]<>",d}"->"s"<>ToString[indno]<>",u}"];

,{indno,1,Length[lst]}];

tmpeval="Do[tracetmp=ST[tracetmp,{"<>StringDrop[tmpstring1,-1]<>"},Sum[GT[tmpt1,{"<>StringDrop[tmpstring2,-1]<>"}],"<>StringDrop[tmpstring4,-1]<>"]],"<>StringDrop[tmpstring3,-1]<>"]";

ToExpression[tmpeval];

];

tracetmp
)

Syntax:   outerT[tensor1,tensor2]

outerT[t1_,t2_]:=(

dimt1=Length[Dimensions[t1[[2]]]];dimt2=Length[Dimensions[t2[[2]]]];
tmpt1=t1;tmpt2=t2;
Clear[outertmp];

MT[outertmp,dimt1+dimt2,Table[d,{ii,1,dimt1+dimt2}]];

tmpstring1="";tmpstring2="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

tmpstring3="";tmpstring4="";

Do[

tmpstring3=tmpstring3<>"{j"<>ToString[indno]<>",d},";
tmpstring4=tmpstring4<>"{j"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt2}];

tmpeval="Do[outertmp=ST[outertmp,{"<>tmpstring1<>StringDrop[tmpstring3,-1]<>"},GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[tmpt2,{"<>StringDrop[tmpstring3,-1]<>"}]],"<>tmpstring2<>StringDrop[tmpstring4,-1]<>"]";

ToExpression[tmpeval];

outertmp

)

scalarDT[t1_]:=(

tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,1,{d}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{{i1,d}},D[tmpt1,ind[[i1]]]],{i1,1,Nind}]";

ToExpression[tmpeval];

covDTtmp

)

covDT[t1_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,dimt1+1,Table[d,{ii,1,dimt1+1}]];

(* Build expression *)

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

Do[
tmpstring4="-Sum[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[Γ,{{REPL,d},{s1,d},{s2,u}}],{s2,1,Nind}]";

tmpstring4=StringReplace[tmpstring4,"i"<>ToString[indno]-> "s2"];

tmpstring4=StringReplace[tmpstring4,"REPL"->"i"<>ToString[indno]];

tmpstring3=tmpstring3<>tmpstring4;

,{indno,1,dimt1}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{"<>tmpstring1<>"{s1,d}},D[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}],ind[[s1]]]"<>tmpstring3<>"],"<>tmpstring2<>"{s1,1,Nind}]";


ToExpression[tmpeval];

covDTtmp

);

.

Buildmetric

Gd=(1/(z^2))(({{-T[v,z], -gvz[v,z], 0, 0}, {-gvz[v,z], gzz[v,z], 0, 0}, {0, 0, S[v,z]^2, 0}, {0, 0, 0, S[v,z]^2}}))

ind={v,z,x1,x2}

Nind=Length[ind]

BuildMetric;

Gd

Gu

detg=Simplify[Det[Gd]]

GRrun2.nb

Clear[Γ]

MT[Γ,3,{d,d,u}];

Clear[R4]

MT[R4,4,{d,d,d,u}];

Clear[R2]

MT[R2,2,{d,d}];

Clear[G2]

MT[G2,2,{d,d}];

Do[Γ=ST[Γ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}},Sum[(1/2)GT[G,{{ind[[i3]],u},{ind[[ss1]],u}}](D[GT[G,{{ind[[i1]],d},{ind[[ss1]],d}}],ind[[i2]]]+D[GT[G,{{ind[[ss1]],d},{ind[[i2]],d}}],ind[[i1]]]-D[GT[G,{{ind[[i1]],d},{ind[[i2]],d}}],ind[[ss1]]]),{ss1,1,Nind}]],{i1,1,Nind},{i2,1,Nind},{i3,1,Nind}]

(* FullSimplify[Γ] *)

Γ=FullSimplify[Γ];

Do[
	R4=ST[R4,
			{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}},
	
			D[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i2]]]-D[GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i1]]]+
			Sum[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i2]],d},{ind[[ss1]],d},{ind[[i4]],u}}]-GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i1]],d},{ind[[ss1]],d},{ind[[i4]],u}}],{ss1,1,Nind}]
		
		]
		,
	{i1,1,Nind},{i2,1,Nind},{i3,1,Nind},{i4,1,Nind}]

(* FullSimplify[R4] *)

Do[R2=ST[R2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			Sum[GT[R4,{{ind[[i1]],d},{ind[[ss1]],u},{ind[[i2]],d},{ind[[ss1]],d}}],{ss1,1,Nind}]
		]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[R2] *)

Ricci=Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}];

Do[G2=ST[G2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			GT[R2,{{ind[[i1]],d},{ind[[i2]],d}}]-(1/2)GT[G,{{ind[[i1]],d},{ind[[i2]],d}}]Ricci]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[G2] *)

(* FullSimplify[Ricci] *)

Einstein equations

R2=Simplify[R2];

Ricci=Simplify[Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}]];

Clear[A];

MT[A,1,{d}];

A=ST[A,{{v,d}},av[v,z]];

A=ST[A,{{z,d}},az[v,z]];

Clear[F];

MT[F,2,{d,d}];

Do[
F=ST[F,{{i1,d},{i2,d}},D[GT[A,{{i1,d}}],ind[[i2]]]-D[GT[A,{{i2,d}}],ind[[i1]]]];
,{i1,1,Nind},{i2,1,Nind}];

trFF=Simplify[trT[outerT[F,F],{{2,4}}]];

Fsqr=Simplify[trT[trFF,{{1,2}}]];

ϕ[v_,z_]:=p[v,z]+I α[v,z]

Clear[Dp];

MT[Dp,1,{d}];

Do[
Dp=ST[Dp,{{i1,d}},D[ϕ[v,z],ind[[i1]]]-I GT[A,{{i1,d}}]ϕ[v,z]];
,{i1,1,Nind}];

ConjDp={{d},Conjugate[Dp[[2]]]};

DpDp=Simplify[ComplexExpand[trT[outerT[ConjDp,Dp],{1,2}]]];

Vterm=Simplify[-2ComplexExpand[Conjugate[ϕ[v,z]]ϕ[v,z]]];

DpDp=(DpDp/.{Re[xx_]->xx,Im[xx_]->0});

sqrtdetg=Simplify[PowerExpand[Sqrt[Det[-Gd]]]];

L=sqrtdetg(Ricci+6-(1/4)Fsqr-DpDp-Vterm);

tmp1=Simplify[Coefficient[L,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]];

tmp2=Simplify[Coefficient[L,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]];

tmp3=Simplify[Coefficient[L,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]];

tmp4=Simplify[Coefficient[L,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]];

GHBT=Simplify[D[tmp1.{Derivative[1,0][gvz][v,z],0,(1/2)Derivative[0,1][gvz][v,z]},v]+D[tmp1.{0,Derivative[0,1][gvz][v,z],(1/2)Derivative[1,0][gvz][v,z]},z]];

GHBT+=Simplify[D[tmp2.{Derivative[1,0][gzz][v,z],0,(1/2)Derivative[0,1][gzz][v,z]},v]+D[tmp2.{0,Derivative[0,1][gzz][v,z],(1/2)Derivative[1,0][gzz][v,z]},z]];

GHBT+=Simplify[D[tmp3.{Derivative[1,0][T][v,z],0,(1/2)Derivative[0,1][T][v,z]},v]+D[tmp3.{0,Derivative[0,1][T][v,z],(1/2)Derivative[1,0][T][v,z]},z]];

GHBT+=Simplify[D[tmp4.{Derivative[1,0][S][v,z],0,(1/2)Derivative[0,1][S][v,z]},v]+D[tmp4.{0,Derivative[0,1][S][v,z],(1/2)Derivative[1,0][S][v,z]},z]];

L2=L-GHBT;

tmp1=Simplify[Coefficient[L2,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]]

tmp2=Simplify[Coefficient[L2,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]]

tmp3=Simplify[Coefficient[L2,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]]

tmp4=Simplify[Coefficient[L2,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]]

EOM[fld_]:=(
tmp=D[L2,fld];
Do[
If[(tmp2=D[fld,ind[[ii]]])=!=0,
tmp-=D[D[L2,tmp2],ind[[ii]]];
]
,{ii,1,Nind}];
tmp
)

eqp=EOM[p[v,z]];

eqα=EOM[α[v,z]];

eqav=EOM[av[v,z]];

eqaz=EOM[az[v,z]];

eqgvz=EOM[gvz[v,z]];

eqgzz=EOM[gzz[v,z]];

eqT=EOM[T[v,z]];

eqS=EOM[S[v,z]];

Clear[gvz,gzz,az]

gvz[v_,z_]:=1

gzz[v_,z_]:=0

az[v_,z_]:=0

eqp2=Simplify[eqp]

eqα2=Simplify[eqα]

eqav2=Simplify[eqav]

eqaz2=Simplify[eqaz]

eqgvz2=Simplify[eqgvz]

eqgzz2=Simplify[eqgzz]

eqT2=Simplify[eqT]

eqS2=Simplify[eqS]

p[v_,z_]:=0

α[v_,z_]:=0

av[v_,z_]:=1-q z

T[v_,z_]:=1-z^3+(1/4)q^2 z^4

S[v_,z_]:=1

Simplify[eqp2]

Simplify[eqα2]

Simplify[eqav2]

Simplify[eqaz2]

Simplify[eqgvz2]

Simplify[eqgzz2]

Simplify[eqT2]

Simplify[eqS2]

Clear[α,p,av,T,S]

p[v_,z_]:=z f[v,z]

α[v_,z_]:=z h[v,z]

av[v_,z_]:=m[v]-q z+z^2a[v,z]

T[v_,z_]:=1+z^2 t[v,z]

S[v_,z_]:=1+z^2 s[v,z]

Evsol=FullSimplify[Solve[{eqgvz2==0,eqS2==0,eqp2==0,eqα2==0,eqaz==0},{Derivative[1,1][f][v,z],Derivative[1,1][h][v,z], Derivative[1,1][a][v,z], Derivative[1,1][s][v,z], Derivative[0,2][t][v,z]}][[1]]]

Momsol=FullSimplify[Solve[{eqT2==0,eqgzz==0,eqav==0},{Derivative[2,0][s][v,z], Derivative[0,2][s][v,z], Derivative[0,2][a][v,z]}][[1]]]

Simplify[D[eqT2,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqav,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqgzz,z]/.D[Evsol,v]/.Evsol/.Momsol]

Const1=Derivative[0,2][s][v,z]-(Derivative[0,2][s][v,z]/.Momsol)

Const2=Derivative[0,2][a][v,z]-(Derivative[0,2][a][v,z]/.Momsol)

Const3=Derivative[2,0][s][v,z]-(Derivative[2,0][s][v,z]/.Momsol)

Eq1=Derivative[1,1][f][v,z]-(Derivative[1,1][f][v,z]/.Evsol)

Eq2=Derivative[1,1][h][v,z]-(Derivative[1,1][h][v,z]/.Evsol)

Eq3=Derivative[1,1][a][v,z]-(Derivative[1,1][a][v,z]/.Evsol)

Eq4=Derivative[1,1][s][v,z]-(Derivative[1,1][s][v,z]/.Evsol)

Eq5=Derivative[0,2][t][v,z]-(Derivative[0,2][t][v,z]/.Evsol)

Clear[deriv];

deriv[X_]:={X[v,z]->X,D[X[v,z],v]->ToExpression[ToString[X]<>"v"],D[X[v,z],z]->ToExpression[ToString[X]<>"z"],D[D[X[v,z],v],z]->ToExpression[ToString[X]<>"vz"],D[X[v,z],{v,2}]->ToExpression[ToString[X]<>"vv"],D[X[v,z],{z,2}]->ToExpression[ToString[X]<>"zz"]};

deriv[j]

repl=Join[{m[v]->m},deriv[f],deriv[a],deriv[h],deriv[t],deriv[s]];

C1str="";

C2str="";

C3str="";

tmp=ToString[CForm[Eq1//.repl//.{Power->power}]];

Eq1str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq2//.repl//.{Power->power}]];

Eq2str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq3//.repl//.{Power->power}]];

Eq3str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq4//.repl//.{Power->power}]];

Eq4str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq5//.repl//.{Power->power}]];

Eq5str=StringReplace[tmp,"*"->".*"]

Export["~/Work/AdSCMTDyn/matlab/eqn.text","\n\neq1 = "<>Eq1str<>" ; \n\neq2 = "<>Eq2str<>" ; \n\neq3 = "<>Eq3str<>" ; \n\neq4 = "<>Eq4str<>" ; \n\neq5 = "<>Eq5str<>" ; \n\nconst1 = "<>C1str<>" ; \n\nconst2 = "<>C2str<>" ; \n\nconst3 = "<>C3str<>" ; \n\n","string"]

Dynamic AdS/CMT - Superconductor - Spectral - Numerial Ansatz II

Print["hmm"]

GRSetup2.nb - modified

(* Latest Version *)

metric={gtt,grr,gθθ,gϕϕ,gzz}

Gd=DiagonalMatrix[metric]

Gd[[1,2]]=grt

Gd[[2,1]]=grt

BuildMetric:=(Gu=FullSimplify[Inverse[Gd]];G={{u,d},IdentityMatrix[Nind]};)

BuildMetric

ind={t,r,θ,ϕ,z}

Nind=Length[ind]

getind[i1_]:=(If[IntegerQ[i1],i1,Position[ind,i1][[1,1]]])

lstind[i1_]:=(lsttmp2={};Do[AppendTo[lsttmp2,getind[i1[[lsttmp1,1]]]],{lsttmp1,1,Length[i1]}];lsttmp2)

up[i1_,i2_]:=(lsttmp2={};Do[AppendTo[lsttmp2,If[i1[[lsttmp1,2]]=!=i2[[lsttmp1]],If[i1[[lsttmp1,2]]===d,d,u],0]],{lsttmp1,1,Length[i1]}];lsttmp2)

MT[t1_,dim_,indpos_]:=(MTtmp1=0;If[Length[indpos]=!=dim,Print["Error in index position"]];Do[MTtmp1=Table[MTtmp1,{MTtmp3,1,Nind}],{MTtmp2,1,dim}];t1={indpos,MTtmp1})

ST[t1_,tind_,val_]:=(STtmp1=t1;If[up[tind,t1[[1]]]=!=Table[0,{Length[tind]}],Print["Error in Index Assignment!!"];];STtmp1[[2]]=ReplacePart[t1[[2]],val,lstind[tind]];STtmp1)

Edited the expression below as 7.0 seems to complain about Sum[x] whereas 6.0 would just return x (ie. if there is no sum index!)

GT[t1_,tind_]:=(GTtmp2=up[tind,t1[[1]]];GTtmp3="";GTtmp4="";GTtmp5="";GTtmp6=t1;
		Do[
			Switch[GTtmp2[[GTtmp1]],
				0,
				(GTtmp3=GTtmp3<>","<>ToString[getind[tind[[GTtmp1,1]]]];),
				u,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gu[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}"),
				d,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gd[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}";)];
			,{GTtmp1,1,Length[t1[[1]]]}];If[GTtmp4===""&&GTtmp5==="",ToExpression["GTtmp6[[2"<>GTtmp3<>"]]"],ToExpression["Sum[GTtmp6[[2"<>GTtmp3<>"]]"<>GTtmp4<>GTtmp5<>"]"]])

Trace, Outer Product and Covariant Derivative! - GREinstein.nb

Syntax:   trT[tensor,{{1,4},{2,3}}]   -  ie. contract 1 and 4 plus 2 and 3 positions

trT[t1_,lst_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[tracetmp];

If[dimt1===2,

tracetmp=Sum[GT[tmpt1,{{s1,d},{s1,u}}],{s1,1,Nind}];

,

MT[tracetmp,dimt1-2Length[lst],Table[d,{ii,1,dimt1-2Length[lst]}]];

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

If[{}===Position[lst,indno],
tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring3=tmpstring3<>"{i"<>ToString[indno]<>",1,Nind},";
];

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",d},";

,{indno,1,dimt1}];

tmpstring4="";

Do[

tmpstring4=tmpstring4<>"{s"<>ToString[indno]<>",1,Nind},";

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,1]]]->"s"<>ToString[indno]];

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,2]]]<>",d}"->"s"<>ToString[indno]<>",u}"];

,{indno,1,Length[lst]}];

tmpeval="Do[tracetmp=ST[tracetmp,{"<>StringDrop[tmpstring1,-1]<>"},Sum[GT[tmpt1,{"<>StringDrop[tmpstring2,-1]<>"}],"<>StringDrop[tmpstring4,-1]<>"]],"<>StringDrop[tmpstring3,-1]<>"]";

ToExpression[tmpeval];

];

tracetmp
)

Syntax:   outerT[tensor1,tensor2]

outerT[t1_,t2_]:=(

dimt1=Length[Dimensions[t1[[2]]]];dimt2=Length[Dimensions[t2[[2]]]];
tmpt1=t1;tmpt2=t2;
Clear[outertmp];

MT[outertmp,dimt1+dimt2,Table[d,{ii,1,dimt1+dimt2}]];

tmpstring1="";tmpstring2="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

tmpstring3="";tmpstring4="";

Do[

tmpstring3=tmpstring3<>"{j"<>ToString[indno]<>",d},";
tmpstring4=tmpstring4<>"{j"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt2}];

tmpeval="Do[outertmp=ST[outertmp,{"<>tmpstring1<>StringDrop[tmpstring3,-1]<>"},GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[tmpt2,{"<>StringDrop[tmpstring3,-1]<>"}]],"<>tmpstring2<>StringDrop[tmpstring4,-1]<>"]";

ToExpression[tmpeval];

outertmp

)

scalarDT[t1_]:=(

tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,1,{d}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{{i1,d}},D[tmpt1,ind[[i1]]]],{i1,1,Nind}]";

ToExpression[tmpeval];

covDTtmp

)

covDT[t1_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,dimt1+1,Table[d,{ii,1,dimt1+1}]];

(* Build expression *)

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

Do[
tmpstring4="-Sum[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[Γ,{{REPL,d},{s1,d},{s2,u}}],{s2,1,Nind}]";

tmpstring4=StringReplace[tmpstring4,"i"<>ToString[indno]-> "s2"];

tmpstring4=StringReplace[tmpstring4,"REPL"->"i"<>ToString[indno]];

tmpstring3=tmpstring3<>tmpstring4;

,{indno,1,dimt1}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{"<>tmpstring1<>"{s1,d}},D[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}],ind[[s1]]]"<>tmpstring3<>"],"<>tmpstring2<>"{s1,1,Nind}]";


ToExpression[tmpeval];

covDTtmp

);

.

Buildmetric

Gd=(1/(z^2))(({{-T[v,z], -gvz[v,z], 0, 0}, {-gvz[v,z], gzz[v,z], 0, 0}, {0, 0, S[v,z]^2, 0}, {0, 0, 0, S[v,z]^2}}))

ind={v,z,x1,x2}

Nind=Length[ind]

BuildMetric;

Gd

Gu

detg=Simplify[Det[Gd]]

GRrun2.nb

Clear[Γ]

MT[Γ,3,{d,d,u}];

Clear[R4]

MT[R4,4,{d,d,d,u}];

Clear[R2]

MT[R2,2,{d,d}];

Clear[G2]

MT[G2,2,{d,d}];

Do[Γ=ST[Γ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}},Sum[(1/2)GT[G,{{ind[[i3]],u},{ind[[ss1]],u}}](D[GT[G,{{ind[[i1]],d},{ind[[ss1]],d}}],ind[[i2]]]+D[GT[G,{{ind[[ss1]],d},{ind[[i2]],d}}],ind[[i1]]]-D[GT[G,{{ind[[i1]],d},{ind[[i2]],d}}],ind[[ss1]]]),{ss1,1,Nind}]],{i1,1,Nind},{i2,1,Nind},{i3,1,Nind}]

(* FullSimplify[Γ] *)

Γ=FullSimplify[Γ];

Do[
	R4=ST[R4,
			{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}},
	
			D[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i2]]]-D[GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i1]]]+
			Sum[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i2]],d},{ind[[ss1]],d},{ind[[i4]],u}}]-GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i1]],d},{ind[[ss1]],d},{ind[[i4]],u}}],{ss1,1,Nind}]
		
		]
		,
	{i1,1,Nind},{i2,1,Nind},{i3,1,Nind},{i4,1,Nind}]

(* FullSimplify[R4] *)

Do[R2=ST[R2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			Sum[GT[R4,{{ind[[i1]],d},{ind[[ss1]],u},{ind[[i2]],d},{ind[[ss1]],d}}],{ss1,1,Nind}]
		]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[R2] *)

Ricci=Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}];

Do[G2=ST[G2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			GT[R2,{{ind[[i1]],d},{ind[[i2]],d}}]-(1/2)GT[G,{{ind[[i1]],d},{ind[[i2]],d}}]Ricci]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[G2] *)

(* FullSimplify[Ricci] *)

Einstein equations

R2=Simplify[R2];

Ricci=Simplify[Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}]];

Clear[A];

MT[A,1,{d}];

A=ST[A,{{v,d}},av[v,z]];

A=ST[A,{{z,d}},az[v,z]];

Clear[F];

MT[F,2,{d,d}];

Do[
F=ST[F,{{i1,d},{i2,d}},D[GT[A,{{i1,d}}],ind[[i2]]]-D[GT[A,{{i2,d}}],ind[[i1]]]];
,{i1,1,Nind},{i2,1,Nind}];

trFF=Simplify[trT[outerT[F,F],{{2,4}}]];

Fsqr=Simplify[trT[trFF,{{1,2}}]];

ϕ[v_,z_]:=p[v,z]+I α[v,z]

Clear[Dp];

MT[Dp,1,{d}];

Do[
Dp=ST[Dp,{{i1,d}},D[ϕ[v,z],ind[[i1]]]-I GT[A,{{i1,d}}]ϕ[v,z]];
,{i1,1,Nind}];

ConjDp={{d},Conjugate[Dp[[2]]]};

DpDp=Simplify[ComplexExpand[trT[outerT[ConjDp,Dp],{1,2}]]];

Vterm=Simplify[-2ComplexExpand[Conjugate[ϕ[v,z]]ϕ[v,z]]];

DpDp=(DpDp/.{Re[xx_]->xx,Im[xx_]->0});

sqrtdetg=Simplify[PowerExpand[Sqrt[Det[-Gd]]]];

L=sqrtdetg(Ricci+6-(1/4)Fsqr-DpDp-Vterm);

tmp1=Simplify[Coefficient[L,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]];

tmp2=Simplify[Coefficient[L,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]];

tmp3=Simplify[Coefficient[L,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]];

tmp4=Simplify[Coefficient[L,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]];

GHBT=Simplify[D[tmp1.{Derivative[1,0][gvz][v,z],0,(1/2)Derivative[0,1][gvz][v,z]},v]+D[tmp1.{0,Derivative[0,1][gvz][v,z],(1/2)Derivative[1,0][gvz][v,z]},z]];

GHBT+=Simplify[D[tmp2.{Derivative[1,0][gzz][v,z],0,(1/2)Derivative[0,1][gzz][v,z]},v]+D[tmp2.{0,Derivative[0,1][gzz][v,z],(1/2)Derivative[1,0][gzz][v,z]},z]];

GHBT+=Simplify[D[tmp3.{Derivative[1,0][T][v,z],0,(1/2)Derivative[0,1][T][v,z]},v]+D[tmp3.{0,Derivative[0,1][T][v,z],(1/2)Derivative[1,0][T][v,z]},z]];

GHBT+=Simplify[D[tmp4.{Derivative[1,0][S][v,z],0,(1/2)Derivative[0,1][S][v,z]},v]+D[tmp4.{0,Derivative[0,1][S][v,z],(1/2)Derivative[1,0][S][v,z]},z]];

L2=L-GHBT;

tmp1=Simplify[Coefficient[L2,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]]

tmp2=Simplify[Coefficient[L2,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]]

tmp3=Simplify[Coefficient[L2,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]]

tmp4=Simplify[Coefficient[L2,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]]

EOM[fld_]:=(
tmp=D[L2,fld];
Do[
If[(tmp2=D[fld,ind[[ii]]])=!=0,
tmp-=D[D[L2,tmp2],ind[[ii]]];
]
,{ii,1,Nind}];
tmp
)

eqp=EOM[p[v,z]];

eqα=EOM[α[v,z]];

eqav=EOM[av[v,z]];

eqaz=EOM[az[v,z]];

eqgvz=EOM[gvz[v,z]];

eqgzz=EOM[gzz[v,z]];

eqT=EOM[T[v,z]];

eqS=EOM[S[v,z]];

Clear[gvz,gzz,az]

gvz[v_,z_]:=1

gzz[v_,z_]:=0

az[v_,z_]:=0

eqp2=Simplify[eqp]

eqα2=Simplify[eqα]

eqav2=Simplify[eqav]

eqaz2=Simplify[eqaz]

eqgvz2=Simplify[eqgvz]

eqgzz2=Simplify[eqgzz]

eqT2=Simplify[eqT]

eqS2=Simplify[eqS]

p[v_,z_]:=0

α[v_,z_]:=0

av[v_,z_]:=1-q z

T[v_,z_]:=1-z^3+(1/4)q^2 z^4

S[v_,z_]:=1

Simplify[eqp2]

Simplify[eqα2]

Simplify[eqav2]

Simplify[eqaz2]

Simplify[eqgvz2]

Simplify[eqgzz2]

Simplify[eqT2]

Simplify[eqS2]

Clear[α,p,av,T,S]

p[v_,z_]:=z f[v,z]

α[v_,z_]:=z h[v,z]

av[v_,z_]:=m[v]+z a[v,z]

T[v_,z_]:=1+z^2 t[v,z]

S[v_,z_]:=1+z^2 s[v,z]

Evsol=FullSimplify[Solve[{eqgvz2==0,eqS2==0,eqp2==0,eqα2==0,eqaz==0},{Derivative[1,1][f][v,z],Derivative[1,1][h][v,z], Derivative[1,1][a][v,z], Derivative[1,1][s][v,z], Derivative[0,2][t][v,z]}][[1]]]

Momsol=FullSimplify[Solve[{eqT2==0,eqgzz==0,eqav==0},{Derivative[2,0][s][v,z], Derivative[0,2][s][v,z], Derivative[0,2][a][v,z]}][[1]]]

Simplify[D[eqT2,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqav,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqgzz,z]/.D[Evsol,v]/.Evsol/.Momsol]

Const1=Derivative[0,2][s][v,z]-(Derivative[0,2][s][v,z]/.Momsol)

Const2=Derivative[0,2][a][v,z]-(Derivative[0,2][a][v,z]/.Momsol)

Const3=Derivative[2,0][s][v,z]-(Derivative[2,0][s][v,z]/.Momsol)

Eq1=Derivative[1,1][f][v,z]-(Derivative[1,1][f][v,z]/.Evsol)

Eq2=Derivative[1,1][h][v,z]-(Derivative[1,1][h][v,z]/.Evsol)

Eq3=Derivative[1,1][a][v,z]-(Derivative[1,1][a][v,z]/.Evsol)

Eq4=Derivative[1,1][s][v,z]-(Derivative[1,1][s][v,z]/.Evsol)

Eq5=Derivative[0,2][t][v,z]-(Derivative[0,2][t][v,z]/.Evsol)

Clear[deriv];

deriv[X_]:={X[v,z]->X,D[X[v,z],v]->ToExpression[ToString[X]<>"v"],D[X[v,z],z]->ToExpression[ToString[X]<>"z"],D[D[X[v,z],v],z]->ToExpression[ToString[X]<>"vz"],D[X[v,z],{v,2}]->ToExpression[ToString[X]<>"vv"],D[X[v,z],{z,2}]->ToExpression[ToString[X]<>"zz"]};

deriv[j]

repl=Join[{m[v]->m},deriv[f],deriv[a],deriv[h],deriv[t],deriv[s]];

tmp=ToString[CForm[Eq1//.repl//.{Power->power}]];

Eq1str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq2//.repl//.{Power->power}]];

Eq2str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq3//.repl//.{Power->power}]];

Eq3str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq4//.repl//.{Power->power}]];

Eq4str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Eq5//.repl//.{Power->power}]];

Eq5str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Const1//.repl//.{Power->power}]];

C1str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Const2//.repl//.{Power->power}]];

C2str=StringReplace[tmp,"*"->".*"]

tmp=ToString[CForm[Const3//.repl//.{Power->power}]];

C3str=StringReplace[tmp,"*"->".*"]

Export["~/Work/AdSCMTDyn2/matlab/eqn.text","\n\neq1 = "<>Eq1str<>" ; \n\neq2 = "<>Eq2str<>" ; \n\neq3 = "<>Eq3str<>" ; \n\neq4 = "<>Eq4str<>" ; \n\neq5 = "<>Eq5str<>" ; \n\nconst1 = "<>C1str<>" ; \n\nconst2 = "<>C2str<>" ; \n\nconst3 = "<>C3str<>" ; \n\n","string"]

f[v_,z_]:=σ[v]+z f2[v]+z^2 f3[v]+z^3 f4[v]+z^4 f5[v];

h[v_,z_]:=z h2[v]+z^2 h3[v]+z^3 h4[v]+z^4 h5[v];

a[v_,z_]:=a1[v]+z a2[v]+z^2 a3[v]+z^3 a4[v]+z^4 a5[v];

t[v_,z_]:=-(1/2) ((σ[v])^2)+z t3[v]+z^2 t4[v]+z^3 t5[v];

s[v_,z_]:=-(1/4) ((σ[v])^2) +z s3[v]+z^2 s4[v]+z^3 s5[v];

tmp1=Simplify[Series[Eq1,{z,0,0}]]

tmp2=Simplify[Series[Eq2,{z,0,0}]]

tmp3=Simplify[Series[Eq3,{z,0,-1}]]

tmp4=Simplify[Series[Eq4,{z,0,-1}]]

tmp5=Simplify[Series[Eq5,{z,0,-1}]]

tmp6=Simplify[Series[Const1,{z,0,-1}]]

tmp7=Simplify[Series[Const2,{z,0,-1}]]

tmp8=Simplify[Series[Const3,{z,0,-1}]]

repl1=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{a1'[v],a2[v],f3[v],h3[v],s3[v]}][[1]]]

repl=Join[ repl1,D[repl1,v]];

tmp1=Simplify[Series[Eq1/.repl,{z,0,1}]]

tmp2=Simplify[Series[Eq2/.repl,{z,0,1}]]

tmp3=Simplify[Series[Eq3/.repl,{z,0,0}]]

tmp4=Simplify[Series[Eq4/.repl,{z,0,0}]]

tmp5=Simplify[Series[Eq5/.repl,{z,0,0}]]

tmp6=Simplify[Series[Const1/.repl,{z,0,0}]]

tmp7=Simplify[Series[Const2/.repl,{z,0,0}]]

tmp8=Simplify[Series[Const3/.repl,{z,0,0}]]

repl2=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{t3'[v],a3[v],t4[v],s4[v],h4[v],f4[v]}][[1]]]

repl=Join[ repl1,D[repl1,v],D[repl1,{v,2}],repl2,D[repl2,v]];

tmp1=Simplify[Series[Eq1//.repl,{z,0,2}]]

tmp2=Simplify[Series[Eq2//.repl,{z,0,2}]]

tmp3=Simplify[Series[Eq3//.repl,{z,0,1}]]

tmp4=Simplify[Series[Eq4//.repl,{z,0,1}]]

tmp5=Simplify[Series[Eq5//.repl,{z,0,1}]]

tmp6=Simplify[Series[Const1//.repl,{z,0,1}]]

tmp7=Simplify[Series[Const2//.repl,{z,0,1}]]

tmp8=Simplify[Series[Const3//.repl,{z,0,1}]]

repl3=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{a4[v],t5[v],s5[v],h5[v],f5[v]}][[1]]]

bceq1=D[a[v,z],v]+2σ[v](D[h[v,z],z]+m[v] σ[v]);

tmp=Simplify[Series[bceq1,{z,0,0}]]

Simplify[tmp/.repl]

bceq2=D[t[v,z],v]+m[v]f[v,z](h[v,z]+z m[v] f[v,z])+f[v,z]D[f[v,z],v]-z D[f[v,z],z]D[f[v,z],v]- z f[v,z] D[f[v,z],{v,2}];

bceq2=D[t[v,z],v]+m[v]σ[v](h[v,z]+z m[v] σ[v])+σ[v]D[f[v,z],v]-z σ[v] σ''[v];

tmp=Simplify[Series[bceq2,{z,0,1}]]

Simplify[tmp/.repl2]

Dynamic AdS/CMT - Superconductor - Spectral - Broken phase

Print["hmm"]

GRSetup2.nb - modified

(* Latest Version *)

metric={gtt,grr,gθθ,gϕϕ,gzz}

Gd=DiagonalMatrix[metric]

Gd[[1,2]]=grt

Gd[[2,1]]=grt

BuildMetric:=(Gu=FullSimplify[Inverse[Gd]];G={{u,d},IdentityMatrix[Nind]};)

BuildMetric

ind={t,r,θ,ϕ,z}

Nind=Length[ind]

getind[i1_]:=(If[IntegerQ[i1],i1,Position[ind,i1][[1,1]]])

lstind[i1_]:=(lsttmp2={};Do[AppendTo[lsttmp2,getind[i1[[lsttmp1,1]]]],{lsttmp1,1,Length[i1]}];lsttmp2)

up[i1_,i2_]:=(lsttmp2={};Do[AppendTo[lsttmp2,If[i1[[lsttmp1,2]]=!=i2[[lsttmp1]],If[i1[[lsttmp1,2]]===d,d,u],0]],{lsttmp1,1,Length[i1]}];lsttmp2)

MT[t1_,dim_,indpos_]:=(MTtmp1=0;If[Length[indpos]=!=dim,Print["Error in index position"]];Do[MTtmp1=Table[MTtmp1,{MTtmp3,1,Nind}],{MTtmp2,1,dim}];t1={indpos,MTtmp1})

ST[t1_,tind_,val_]:=(STtmp1=t1;If[up[tind,t1[[1]]]=!=Table[0,{Length[tind]}],Print["Error in Index Assignment!!"];];STtmp1[[2]]=ReplacePart[t1[[2]],val,lstind[tind]];STtmp1)

Edited the expression below as 7.0 seems to complain about Sum[x] whereas 6.0 would just return x (ie. if there is no sum index!)

GT[t1_,tind_]:=(GTtmp2=up[tind,t1[[1]]];GTtmp3="";GTtmp4="";GTtmp5="";GTtmp6=t1;
		Do[
			Switch[GTtmp2[[GTtmp1]],
				0,
				(GTtmp3=GTtmp3<>","<>ToString[getind[tind[[GTtmp1,1]]]];),
				u,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gu[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}"),
				d,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gd[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}";)];
			,{GTtmp1,1,Length[t1[[1]]]}];If[GTtmp4===""&&GTtmp5==="",ToExpression["GTtmp6[[2"<>GTtmp3<>"]]"],ToExpression["Sum[GTtmp6[[2"<>GTtmp3<>"]]"<>GTtmp4<>GTtmp5<>"]"]])

Trace, Outer Product and Covariant Derivative! - GREinstein.nb

Syntax:   trT[tensor,{{1,4},{2,3}}]   -  ie. contract 1 and 4 plus 2 and 3 positions

trT[t1_,lst_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[tracetmp];

If[dimt1===2,

tracetmp=Sum[GT[tmpt1,{{s1,d},{s1,u}}],{s1,1,Nind}];

,

MT[tracetmp,dimt1-2Length[lst],Table[d,{ii,1,dimt1-2Length[lst]}]];

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

If[{}===Position[lst,indno],
tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring3=tmpstring3<>"{i"<>ToString[indno]<>",1,Nind},";
];

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",d},";

,{indno,1,dimt1}];

tmpstring4="";

Do[

tmpstring4=tmpstring4<>"{s"<>ToString[indno]<>",1,Nind},";

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,1]]]->"s"<>ToString[indno]];

tmpstring2=StringReplace[tmpstring2,"i"<>ToString[lst[[indno,2]]]<>",d}"->"s"<>ToString[indno]<>",u}"];

,{indno,1,Length[lst]}];

tmpeval="Do[tracetmp=ST[tracetmp,{"<>StringDrop[tmpstring1,-1]<>"},Sum[GT[tmpt1,{"<>StringDrop[tmpstring2,-1]<>"}],"<>StringDrop[tmpstring4,-1]<>"]],"<>StringDrop[tmpstring3,-1]<>"]";

ToExpression[tmpeval];

];

tracetmp
)

Syntax:   outerT[tensor1,tensor2]

outerT[t1_,t2_]:=(

dimt1=Length[Dimensions[t1[[2]]]];dimt2=Length[Dimensions[t2[[2]]]];
tmpt1=t1;tmpt2=t2;
Clear[outertmp];

MT[outertmp,dimt1+dimt2,Table[d,{ii,1,dimt1+dimt2}]];

tmpstring1="";tmpstring2="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";
tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

tmpstring3="";tmpstring4="";

Do[

tmpstring3=tmpstring3<>"{j"<>ToString[indno]<>",d},";
tmpstring4=tmpstring4<>"{j"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt2}];

tmpeval="Do[outertmp=ST[outertmp,{"<>tmpstring1<>StringDrop[tmpstring3,-1]<>"},GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[tmpt2,{"<>StringDrop[tmpstring3,-1]<>"}]],"<>tmpstring2<>StringDrop[tmpstring4,-1]<>"]";

ToExpression[tmpeval];

outertmp

)

scalarDT[t1_]:=(

tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,1,{d}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{{i1,d}},D[tmpt1,ind[[i1]]]],{i1,1,Nind}]";

ToExpression[tmpeval];

covDTtmp

)

covDT[t1_]:=(

dimt1=Length[Dimensions[t1[[2]]]];
tmpt1=t1;
Clear[covDTtmp];
MT[covDTtmp,dimt1+1,Table[d,{ii,1,dimt1+1}]];

(* Build expression *)

tmpstring1="";tmpstring2="";tmpstring3="";

Do[

tmpstring1=tmpstring1<>"{i"<>ToString[indno]<>",d},";

tmpstring2=tmpstring2<>"{i"<>ToString[indno]<>",1,Nind},";

,{indno,1,dimt1}];

Do[
tmpstring4="-Sum[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}]GT[Γ,{{REPL,d},{s1,d},{s2,u}}],{s2,1,Nind}]";

tmpstring4=StringReplace[tmpstring4,"i"<>ToString[indno]-> "s2"];

tmpstring4=StringReplace[tmpstring4,"REPL"->"i"<>ToString[indno]];

tmpstring3=tmpstring3<>tmpstring4;

,{indno,1,dimt1}];

tmpeval="Do[covDTtmp=ST[covDTtmp,{"<>tmpstring1<>"{s1,d}},D[GT[tmpt1,{"<>StringDrop[tmpstring1,-1]<>"}],ind[[s1]]]"<>tmpstring3<>"],"<>tmpstring2<>"{s1,1,Nind}]";


ToExpression[tmpeval];

covDTtmp

);

.

Buildmetric

Gd=(1/(z^2))(({{-T[v,z], -gvz[v,z], 0, 0}, {-gvz[v,z], gzz[v,z], 0, 0}, {0, 0, S[v,z]^2, 0}, {0, 0, 0, S[v,z]^2}}))

ind={v,z,x1,x2}

Nind=Length[ind]

BuildMetric;

Gd

Gu

detg=Simplify[Det[Gd]]

GRrun2.nb

Clear[Γ]

MT[Γ,3,{d,d,u}];

Clear[R4]

MT[R4,4,{d,d,d,u}];

Clear[R2]

MT[R2,2,{d,d}];

Clear[G2]

MT[G2,2,{d,d}];

Do[Γ=ST[Γ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}},Sum[(1/2)GT[G,{{ind[[i3]],u},{ind[[ss1]],u}}](D[GT[G,{{ind[[i1]],d},{ind[[ss1]],d}}],ind[[i2]]]+D[GT[G,{{ind[[ss1]],d},{ind[[i2]],d}}],ind[[i1]]]-D[GT[G,{{ind[[i1]],d},{ind[[i2]],d}}],ind[[ss1]]]),{ss1,1,Nind}]],{i1,1,Nind},{i2,1,Nind},{i3,1,Nind}]

(* FullSimplify[Γ] *)

Γ=FullSimplify[Γ];

Do[
	R4=ST[R4,
			{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}},
	
			D[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i2]]]-D[GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i1]]]+
			Sum[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i2]],d},{ind[[ss1]],d},{ind[[i4]],u}}]-GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i1]],d},{ind[[ss1]],d},{ind[[i4]],u}}],{ss1,1,Nind}]
		
		]
		,
	{i1,1,Nind},{i2,1,Nind},{i3,1,Nind},{i4,1,Nind}]

(* FullSimplify[R4] *)

Do[R2=ST[R2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			Sum[GT[R4,{{ind[[i1]],d},{ind[[ss1]],u},{ind[[i2]],d},{ind[[ss1]],d}}],{ss1,1,Nind}]
		]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[R2] *)

Ricci=Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}];

Do[G2=ST[G2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			GT[R2,{{ind[[i1]],d},{ind[[i2]],d}}]-(1/2)GT[G,{{ind[[i1]],d},{ind[[i2]],d}}]Ricci]
	,{i1,1,Nind},{i2,1,Nind}]

(* FullSimplify[G2] *)

(* FullSimplify[Ricci] *)

Einstein equations

R2=Simplify[R2];

Ricci=Simplify[Sum[GT[R2,{{ind[[ss1]],u},{ind[[ss1]],d}}],{ss1,1,Nind}]];

Clear[A];

MT[A,1,{d}];

A=ST[A,{{v,d}},av[v,z]];

A=ST[A,{{z,d}},az[v,z]];

Clear[F];

MT[F,2,{d,d}];

Do[
F=ST[F,{{i1,d},{i2,d}},D[GT[A,{{i1,d}}],ind[[i2]]]-D[GT[A,{{i2,d}}],ind[[i1]]]];
,{i1,1,Nind},{i2,1,Nind}];

trFF=Simplify[trT[outerT[F,F],{{2,4}}]];

Fsqr=Simplify[trT[trFF,{{1,2}}]];

ϕ[v_,z_]:=p[v,z]+I α[v,z]

Clear[Dp];

MT[Dp,1,{d}];

Do[
Dp=ST[Dp,{{i1,d}},D[ϕ[v,z],ind[[i1]]]-I GT[A,{{i1,d}}]ϕ[v,z]];
,{i1,1,Nind}];

ConjDp={{d},Conjugate[Dp[[2]]]};

DpDp=Simplify[ComplexExpand[trT[outerT[ConjDp,Dp],{1,2}]]];

Vterm=Simplify[-2ComplexExpand[Conjugate[ϕ[v,z]]ϕ[v,z]]];

DpDp=(DpDp/.{Re[xx_]->xx,Im[xx_]->0});

sqrtdetg=Simplify[PowerExpand[Sqrt[Det[-Gd]]]];

L=sqrtdetg(Ricci+6-(1/4)Fsqr-DpDp-Vterm)

tmp1=Simplify[Coefficient[L,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]];

tmp2=Simplify[Coefficient[L,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]];

tmp3=Simplify[Coefficient[L,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]];

tmp4=Simplify[Coefficient[L,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]];

GHBT=Simplify[D[tmp1.{Derivative[1,0][gvz][v,z],0,(1/2)Derivative[0,1][gvz][v,z]},v]+D[tmp1.{0,Derivative[0,1][gvz][v,z],(1/2)Derivative[1,0][gvz][v,z]},z]];

GHBT+=Simplify[D[tmp2.{Derivative[1,0][gzz][v,z],0,(1/2)Derivative[0,1][gzz][v,z]},v]+D[tmp2.{0,Derivative[0,1][gzz][v,z],(1/2)Derivative[1,0][gzz][v,z]},z]];

GHBT+=Simplify[D[tmp3.{Derivative[1,0][T][v,z],0,(1/2)Derivative[0,1][T][v,z]},v]+D[tmp3.{0,Derivative[0,1][T][v,z],(1/2)Derivative[1,0][T][v,z]},z]];

GHBT+=Simplify[D[tmp4.{Derivative[1,0][S][v,z],0,(1/2)Derivative[0,1][S][v,z]},v]+D[tmp4.{0,Derivative[0,1][S][v,z],(1/2)Derivative[1,0][S][v,z]},z]];

L2=L-GHBT;

tmp1=Simplify[Coefficient[L2,{Derivative[2,0][gvz][v,z],Derivative[0,2][gvz][v,z],Derivative[1,1][gvz][v,z]}]]

tmp2=Simplify[Coefficient[L2,{Derivative[2,0][gzz][v,z],Derivative[0,2][gzz][v,z],Derivative[1,1][gzz][v,z]}]]

tmp3=Simplify[Coefficient[L2,{Derivative[2,0][T][v,z],Derivative[0,2][T][v,z],Derivative[1,1][T][v,z]}]]

tmp4=Simplify[Coefficient[L2,{Derivative[2,0][S][v,z],Derivative[0,2][S][v,z],Derivative[1,1][S][v,z]}]]

EOM[fld_]:=(
tmp=D[L2,fld];
Do[
If[(tmp2=D[fld,ind[[ii]]])=!=0,
tmp-=D[D[L2,tmp2],ind[[ii]]];
]
,{ii,1,Nind}];
tmp
)

eqp=EOM[p[v,z]];

eqα=EOM[α[v,z]];

eqav=EOM[av[v,z]];

eqaz=EOM[az[v,z]];

eqgvz=EOM[gvz[v,z]];

eqgzz=EOM[gzz[v,z]];

eqT=EOM[T[v,z]];

eqS=EOM[S[v,z]];

Clear[gvz,gzz,az]

gvz[v_,z_]:=1

gzz[v_,z_]:=0

az[v_,z_]:=0

eqp2=Simplify[eqp]

eqα2=Simplify[eqα]

eqav2=Simplify[eqav]

eqaz2=Simplify[eqaz]

eqgvz2=Simplify[eqgvz]

eqgzz2=Simplify[eqgzz]

eqT2=Simplify[eqT]

eqS2=Simplify[eqS]

Evsol=Simplify[Solve[{eqgvz2==0,eqS2==0,eqp2==0,eqα2==0,eqaz==0},{Derivative[1,1][p][v,z],Derivative[1,1][α][v,z], Derivative[1,1][av][v,z], Derivative[1,1][S][v,z], Derivative[0,2][T][v,z]}][[1]]]

Momsol=Simplify[Solve[{eqT2==0,eqgzz==0,eqav==0},{Derivative[2,0][S][v,z], Derivative[0,2][S][v,z], Derivative[0,2][av][v,z]}][[1]]]

Simplify[D[eqT2,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqav,v]/.D[Evsol,z]/.Evsol/.Momsol]

Simplify[D[eqgzz,z]/.D[Evsol,v]/.Evsol/.Momsol]

av[v_,z_]:=a[v,z]

α[v_,z_]:=h[v,z]

Const1=Derivative[0,2][S][v,z]-(Derivative[0,2][S][v,z]/.Momsol)

Const2=Derivative[0,2][a][v,z]-(Derivative[0,2][a][v,z]/.Momsol)

Const3=Derivative[2,0][S][v,z]-(Derivative[2,0][S][v,z]/.Momsol)

Eq1=Derivative[1,1][p][v,z]-(Derivative[1,1][p][v,z]/.Evsol)

Eq2=Derivative[1,1][h][v,z]-(Derivative[1,1][h][v,z]/.Evsol)

Eq3=Derivative[1,1][a][v,z]-(Derivative[1,1][a][v,z]/.Evsol)

Eq4=Derivative[1,1][S][v,z]-(Derivative[1,1][S][v,z]/.Evsol)

Eq5=Derivative[0,2][T][v,z]-(Derivative[0,2][T][v,z]/.Evsol)

Clear[p,h,a,T,S]

pos = 1;

p[v_,z_]:=(z p1[v]+z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 a2[v]+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1+z^2T2[v]+z^3 T3[v]+z^4 T4[v]+z^5 T5[v];

S[v_,z_]:=1+z^2S2[v]+z^3 S3[v]+z^4 S4[v]+z^5 S5[v];

tmp1=Series[Const1,{z,0,pos-1}]

tmp2=Series[Const2,{z,0,pos-1}]

tmp3=Series[Const3,{z,0,pos-1}]

tmp4=Series[Eq1,{z,0,pos-1}]

tmp5=Series[Eq2,{z,0,pos-1}]

tmp6=Series[Eq3,{z,0,pos-1}]

tmp7=Series[Eq4,{z,0,pos-1}]

tmp8=Series[Eq5,{z,0,pos-1}]

repl=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{a1'[v],S2[v],T2[v],a2[v]}][[1]]]

pos = 1;

p[v_,z_]:=(z p1[v]+z^2 p2[v]+z^3 p3[v]+z^4 p4[v]+z^5 p5[v]);

h[v_,z_]:=(z^2 h2[v]+z^3 h3[v]+z^4 h4[v]+z^5 h5[v]);

a[v_,z_]:=a0[v]+z a1[v]+z^2 (-h2[v] p1[v])+z^3 a3[v]+z^4 a4[v]+z^5 a5[v];

T[v_,z_]:=1+z^2(-(1/2) ((p1[v])^2))+z^3 T3[v]+z^4 T4[v]+z^5 T5[v];

S[v_,z_]:=1+z^2(-(1/4) ((p1[v])^2))+z^3 S3[v]+z^4 S4[v]+z^5 S5[v];

pos=2;

repltmp={a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])};

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl=Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{h3[v],p3[v],S3[v],a3[v]}][[1]]

pos=3;

repltmp=Join[{a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])},repl,D[repl,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl2=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T3'[v],a4[v],h4[v],p4[v],S4[v],T4[v]}][[1]]]

pos=4;

repltmp=Join[{a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])},repl,D[repl,v],D[repl,{v,2}],repl2,D[repl2,v]];

tmp1=Simplify[Series[Const1//.repltmp,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2//.repltmp,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3//.repltmp,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1//.repltmp,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2//.repltmp,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3//.repltmp,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4//.repltmp,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5//.repltmp,{z,0,pos-1}]]

repl3=Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}==0,{T5[v],a5[v],h5[v],p5[v],S5[v],T5[v]}][[1]]]

repltmp=Join[{a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v])},repl,D[repl,v],D[repl,{v,2}],D[repl,{v,3}],repl2,D[repl2,v],D[repl2,{v,2}],repl3,D[repl3,v]];

Simplify[Series[p[v,z],{z,0,3}]/.repltmp]

Simplify[Series[h[v,z],{z,0,3}]/.repltmp]

Simplify[Series[a[v,z],{z,0,4}]/.repltmp]

Simplify[Series[T[v,z],{z,0,4}]/.repltmp]

Simplify[Series[S[v,z],{z,0,4}]/.repltmp]

Normal[Simplify[Series[p[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[h[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[a[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[T[v,z],{z,0,5}]/.repltmp]]

Normal[Simplify[Series[S[v,z],{z,0,5}]/.repltmp]]

p[v_,z_]:=z p1[v]+(z^2) p2[v]+(z^3) (a0[v] h2[v]+(((p1[v])^3)/2)+p2'[v])+(1/6) (z^4) (a1[v] (3 h2[v]-2 a0[v] p1[v])-4 ((a0[v])^2) p2[v]+5 ((p1[v])^2) p2[v]-p1[v] T3[v]+4 h2[v] a0'[v]+8 a0[v] h2'[v]+4 ((p1[v])^2) p1'[v]+4 p2''[v])+(1/12) (z^5) (-4 ((a0[v])^3) h2[v]-3 ((a1[v])^2) p1[v]-8 ((h2[v])^2) p1[v]+4 ((p1[v])^5)+6 p1[v] ((p2[v])^2)-4 p2[v] T3[v]+12 a0'[v] h2'[v]+a1[v] (-2 p1[v] a0'[v]+8 h2'[v])+4 p1[v] p2[v] p1'[v]-T3[v] p1'[v]+8 p1[v] ((p1'[v])^2)+4 ((a0[v])^2) (2 ((p1[v])^3)-3 p2'[v])+13 ((p1[v])^2) p2'[v]+4 h2[v] a0''[v]+a0[v] (17 h2[v] ((p1[v])^2)-4 (3 p2[v] a0'[v]+a1[v] (2 p2[v]+p1'[v])-3 h2''[v]))+3 ((p1[v])^2) p1''[v]+4 Derivative[3][p2][v]);

h[v_,z_]:=(z^2) h2[v]+(z^3) (-(1/2) a1[v] p1[v]-a0[v] p2[v]+h2'[v])+(1/6) (z^4) (-4 ((a0[v])^2) h2[v]+9 h2[v] ((p1[v])^2)-3 a1[v] p2[v]-4 p2[v] a0'[v]-2 a1[v] p1'[v]+a0[v] (3 ((p1[v])^3)-8 p2'[v])+4 h2''[v])+(1/12) (z^5) (4 ((a0[v])^3) p2[v]+14 h2[v] p1[v] p2[v]-4 h2[v] T3[v]+3 ((p1[v])^3) a0'[v]+2 ((a0[v])^2) (a1[v] p1[v]-6 h2'[v])+16 ((p1[v])^2) h2'[v]+20 h2[v] p1[v] p1'[v]-12 a0'[v] p2'[v]-4 p2[v] a0''[v]-2 a1[v] (3 ((p1[v])^3)+4 p2'[v]+p1''[v])+a0[v] (-8 a1[v] h2[v]+p1[v] T3[v]+((p1[v])^2) (-4 p2[v]+9 p1'[v])-12 (h2[v] a0'[v]+p2''[v]))+4 Derivative[3][h2][v]);

a[v_,z_]:=a0[v]+z a1[v]-(z^2) h2[v] p1[v]+(1/6) (z^3) p1[v] (3 a1[v] p1[v]+4 a0[v] p2[v]-4 h2'[v])+(1/12) (z^4) (4 ((a0[v])^2) h2[v] p1[v]+a0[v] (2 ((h2[v])^2)-3 ((p1[v])^4)+2 ((p2[v])^2)+8 p1[v] p2'[v])+2 (2 p1[v] p2[v] a0'[v]-p2[v] h2'[v]+a1[v] p1[v] (3 p2[v]+p1'[v])+h2[v] (-5 ((p1[v])^3)+p2'[v])-2 p1[v] h2''[v]))+(1/30) (z^5) (-4 ((a0[v])^3) p1[v] p2[v]-24 h2[v] ((p1[v])^2) p2[v]+3 h2[v] p1[v] T3[v]+4 ((h2[v])^2) a0'[v]-3 ((p1[v])^4) a0'[v]+4 ((p2[v])^2) a0'[v]+12 ((a0[v])^2) p1[v] h2'[v]-19 ((p1[v])^3) h2'[v]-16 h2[v] ((p1[v])^2) p1'[v]+12 p1[v] a0'[v] p2'[v]+4 p1[v] p2[v] a0''[v]-4 p2[v] h2''[v]+a1[v] (5 ((h2[v])^2)+9 a0[v] h2[v] p1[v]-2 ((a0[v])^2) ((p1[v])^2)+10 ((p1[v])^4)+5 ((p2[v])^2)+2 p2[v] p1'[v]+11 p1[v] p2'[v]+2 p1[v] p1''[v])+4 h2[v] p2''[v]+a0[v] (-((p1[v])^2) T3[v]+((p1[v])^3) (4 p2[v]-9 p1'[v])+8 (h2[v] h2'[v]+p2[v] p2'[v])+12 p1[v] (h2[v] a0'[v]+p2''[v]))-4 p1[v] Derivative[3][h2][v]);

T[v_,z_]:=1-(1/2) (z^2) ((p1[v])^2)+(z^3) T3[v]+(1/12) (z^4) (3 ((a1[v])^2)-2 (2 ((h2[v])^2)+6 a0[v] h2[v] p1[v]+((p1[v])^4)+2 p2[v] (p2[v]-2 p1'[v])+2 p1[v] p2'[v]))+(1/6) (z^5) (a1[v] p1[v] (-3 h2[v]+2 a0[v] p1[v])+4 ((a0[v])^2) p1[v] p2[v]-2 ((p1[v])^3) p2[v]+((p1[v])^2) T3[v]-h2[v] p1[v] a0'[v]-2 h2[v] h2'[v]+a0[v] (-5 p1[v] h2'[v]+3 h2[v] p1'[v])-2 p2[v] p2'[v]+3 p1'[v] p2'[v]-p1[v] p2''[v]);

S[v_,z_]:=1-(1/4) (z^2) ((p1[v])^2)-(1/3) (z^3) p1[v] p2[v]+(1/96) (z^4) (-16 ((h2[v])^2)-24 a0[v] h2[v] p1[v]-11 ((p1[v])^4)-16 ((p2[v])^2)-24 p1[v] p2'[v])+(1/60) (z^5) (a1[v] p1[v] (3 h2[v]+4 a0[v] p1[v])+8 ((a0[v])^2) p1[v] p2[v]-17 ((p1[v])^3) p2[v]+2 ((p1[v])^2) T3[v]-8 h2[v] p1[v] a0'[v]-18 h2[v] h2'[v]-16 a0[v] p1[v] h2'[v]-8 ((p1[v])^3) p1'[v]-18 p2[v] p2'[v]-8 p1[v] p2''[v]);

pos=4;

tmp1=Simplify[Series[Const1,{z,0,pos-1}]]

tmp2=Simplify[Series[Const2,{z,0,pos-1}]]

tmp3=Simplify[Series[Const3,{z,0,pos-1}]]

tmp4=Simplify[Series[Eq1,{z,0,pos-1}]]

tmp5=Simplify[Series[Eq2,{z,0,pos-1}]]

tmp6=Simplify[Series[Eq3,{z,0,pos-1}]]

tmp7=Simplify[Series[Eq4,{z,0,pos-1}]]

tmp8=Simplify[Series[Eq5,{z,0,pos-1}]]

crepl={a1'[v]->-2 p1[v] (h2[v]+a0[v] p1[v]),T3'[v]->-p1[v] (a0[v] h2[v]+((a0[v])^2) p1[v]+p2'[v]-p1''[v])};

Simplify[tmp1]

Simplify[tmp2]

Simplify[tmp3/.crepl]

Simplify[tmp4/.crepl]

Simplify[tmp5/.crepl]

Simplify[tmp6/.crepl]

Simplify[tmp7]

Simplify[tmp8]

Clear[p,h,a,T,S]

p[v_,z_]:=p[z]

h[v_,z_]:=h[z]

a[v_,z_]:=a[z]

T[v_,z_]:=T[z]

S[v_,z_]:=S[z]

ell1=Simplify[Eq1]

ell2=Simplify[Eq2]

ellc1=Simplify[Eq3]

ellc2=Simplify[Eq4]

ell3=Simplify[Eq5]

ell4=Simplify[Const1]

ell5=Simplify[Const2]

ellc3=Simplify[Const3]

ell1=a[z] (h'[z]+h[z] (-(1/z)+((S'[z])/(S[z]))))+(1/2) (-((2 p[z])/(z^2))+h[z] a'[z]+((2 T[z] p'[z])/z)-((2 T[z] p'[z] S'[z])/(S[z]))-p'[z] T'[z]-T[z] p''[z]);

ell2=-(1/(2 (z^2) S[z]))(2 (z^2) T[z] h'[z] S'[z]+2 z a[z] (z S[z] p'[z]+p[z] (-S[z]+z S'[z]))+S[z] (2 h[z]+z (z p[z] a'[z]-2 T[z] h'[z]+z h'[z] T'[z]+z T[z] h''[z])));

ellc1=((2 (a[z] (((h[z])^2)+((p[z])^2))+T[z] (p[z] h'[z]-h[z] p'[z])))/(z^2));

ellc2=-((4 (z^2) T[z] ((S'[z])^2)+4 z S[z] S'[z] (-4 T[z]+z T'[z])-((S[z])^2) (12+4 ((h[z])^2)+4 ((p[z])^2)-(z^4) ((a'[z])^2)+2 T[z] (-6+(z^2) ((h'[z])^2)+(z^2) ((p'[z])^2))+4 z T'[z]))/(8 (z^2) S[z]));

ell3=-((2 z T[z] ((S'[z])^2)+((S[z])^2) ((z^3) ((a'[z])^2)+2 (z a[z] (-p[z] h'[z]+h[z] p'[z])-z T[z] (((h'[z])^2)+((p'[z])^2))+T'[z]))-2 S[z] T[z] (2 S'[z]+z S''[z]))/(z ((S[z])^2)))+T''[z];

ell4=(1/2) S[z] (((h'[z])^2)+((p'[z])^2))+S''[z];

ell5=((2 p[z] h'[z])/(z^2))-((2 h[z] p'[z])/(z^2))+((2 a'[z] S'[z])/(S[z]))+a''[z];

ellc3=(1/(8 (z^2) S[z]))(4 (z^2) ((a[z])^2) (((h[z])^2)+((p[z])^2)) ((S[z])^2)+4 (z^2) a[z] ((S[z])^2) T[z] (p[z] h'[z]-h[z] p'[z])+T[z] (-4 (z^2) T[z] ((S'[z])^2)-4 z S[z] S'[z] (-4 T[z]+z T'[z])+((S[z])^2) (12+4 ((h[z])^2)+4 ((p[z])^2)-12 T[z]-(z^4) ((a'[z])^2)+2 (z^2) T[z] ((h'[z])^2)+2 (z^2) T[z] ((p'[z])^2)+4 z T'[z])));

sol=FullSimplify[Solve[{ell1,ell2,ell3,ell4,ell5}==0,{p''[z],h''[z],a''[z],T''[z],S''[z]}][[1]]]

Expand[ellc1]

Expand[ellc2]

Expand[ellc3]

sol2=FullSimplify[Solve[{ellc1,ellc2}==0,{T'[z],p'[z]}][[1]]]

Simplify[ellc3/.sol2]

Simplify[(D[ellc1,z]+2((S'[z])/(S[z])) ellc1)/.sol]

Simplify[D[S[z]^2ellc1,z]/.sol]

Simplify[(D[ellc2,z]-((S'[z])/(S[z])) ellc2+(2/z) ellc2)/.sol]

Simplify[D[((z^2)/(S[z]))ellc2,z]/.sol]

Print["\n\nSo, if c1 and c2 are zero, and they obey simple decoupled eqns, then c3 also zero!\n\n"]

p[z_]:=(z^2 p2+z^3 p3+z^4 p4+z^5 p5+z^6 p6+z^7 p7+z^8 p8+z^9 p9);

h[z_]:=(z^2 h2+z^3 h3+z^4 h4+z^5 h5+z^6 h6+z^7 h7+z^8 h8+z^9 h9);

a[z_]:=a0+z a1+z^2 a2+z^3 a3+z^4 a4+z^5 a5+z^6 a6+z^7 a7+z^8 a8+z^9 a9;

T[z_]:=1+z^2T2+z^3 T3+z^4 T4+z^5 T5+z^6 T6+z^7 T7+z^8 T8+z^9 T9;

S[z_]:=1+z^2S2+z^3 S3+z^4 S4+z^5 S5+z^6 S6+z^7 S7+z^8 S8+z^9 S9;

runeq:=(

tmp1=Series[ellc1,{z,0,pos-1}];

tmp2=Series[ellc2,{z,0,pos-1}];

tmp3=Series[ellc3,{z,0,pos-1}];

tmp4=Series[ell1,{z,0,pos-1}];

tmp5=Series[ell2,{z,0,pos-1}];

tmp6=Series[ell3,{z,0,pos-1}];

tmp7=Series[ell4,{z,0,pos-1}];

tmp8=Series[ell5,{z,0,pos-1}];

)

pos=1;

runeq

tmplist={tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8}

Solve[tmplist==0,{T2,S2,a2}]

tmp4=Series[ell1,{z,0,pos-1}]

tmp5=Series[ell2,{z,0,pos-1}]

tmp6=Series[ell3,{z,0,pos-1}]

tmp7=Series[ell4,{z,0,pos-1}]

tmp8=Series[ell5,{z,0,pos-1}]

pos = 2;

p[z_]:=(z^2 p2+z^3 p3+z^4 p4+z^5 p5);

h[z_]:=(z^2 h2+z^3 h3+z^4 h4+z^5 h5);

a[z_]:=a0+z a1+z^3 a3+z^4 a4+z^5 a5;

T[z_]:=1+z^3 T3+z^4 T4+z^5 T5;

S[z_]:=1+z^3 S3+z^4 S4+z^5 S5;

tmp1=Series[ellc1,{z,0,pos-1}]

tmp2=Series[ellc2,{z,0,pos-1}]

tmp3=Series[ellc3,{z,0,pos-1}]

tmp4=Series[ell1,{z,0,pos-1}]

tmp5=Series[ell2,{z,0,pos-1}]

tmp6=Series[ell3,{z,0,pos-1}]

tmp7=Series[ell4,{z,0,pos-1}]

tmp8=Series[ell5,{z,0,pos-1}]

pos = 3;

p[z_]:=(z^2 p2+z^3 a0 h2+z^4 p4+z^5 p5);

h[z_]:=(z^2 h2-z^3 a0 p2+z^4 h4+z^5 h5);

a[z_]:=a0+z a1+z^4 a4+z^5 a5;

T[z_]:=1+z^3 T3+z^4 T4+z^5 T5;

S[z_]:=1+z^4 S4+z^5 S5;

tmp1=Series[ellc1,{z,0,pos-1}]

tmp2=Series[ellc2,{z,0,pos-1}]

tmp3=Series[ellc3,{z,0,pos-1}]

tmp4=Series[ell1,{z,0,pos-1}]

tmp5=Series[ell2,{z,0,pos-1}]

tmp6=Series[ell3,{z,0,pos-1}]

tmp7=Series[ell4,{z,0,pos-1}]

tmp8=Series[ell5,{z,0,pos-1}]

pos = 4;

p[z_]:=(z^2 p2+z^3 a0 h2+z^4 (((a1 h2)/2)-(2/3) (a0^2) p2)+z^5 p5);

h[z_]:=(z^2 h2-z^3 a0 p2-z^4 (((a1 p2)/2)+(2/3) (a0^2) h2)+z^5 h5);

a[z_]:=a0+z a1+(1/6)z^4 a0( (h2^2)+ (p2^2))+z^5 a5;

T[z_]:=1+z^3 T3+z^4 (((a1^2)/4)-((((h2^2)+(p2^2)))/3))+z^5 T5;

S[z_]:=1-(1/6)z^4 ((h2^2)+(p2^2))+z^5 S5;

tmp1=Series[ellc1,{z,0,pos-1}]

tmp2=Series[ellc2,{z,0,pos-1}]

tmp3=Series[ellc3,{z,0,pos-1}]

tmp4=Series[ell1,{z,0,pos-1}]

tmp5=Series[ell2,{z,0,pos-1}]

tmp6=Series[ell3,{z,0,pos-1}]

tmp7=Series[ell4,{z,0,pos-1}]

tmp8=Series[ell5,{z,0,pos-1}]

pos = 4;

p[z_]:=(z^2 p2+z^3 a0 h2+z^4 (((a1 h2)/2)-(2/3) (a0^2) p2)-z^5 ((1/3) (a0^3) h2+(2/3) a0 a1 p2+(1/3) p2 T3));

h[z_]:=(z^2 h2-z^3 a0 p2-z^4 (((a1 p2)/2)+(2/3) (a0^2) h2)+z^5 ((1/3) (a0^3) p2-(2/3) a0 a1 h2-(1/3) h2 T3));

a[z_]:=a0+z a1+(1/6)z^4 a0( (h2^2)+ (p2^2))+(1/6) z^5 a1 ( (h2^2)+ (p2^2));

T[z_]:=1+z^3 T3+z^4 (((a1^2)/4)-((((h2^2)+(p2^2)))/3))+z^5 *0;

S[z_]:=1-(1/6)z^4 ((h2^2)+(p2^2))+z^5 *0;

tmp1=Series[ellc1,{z,0,pos-1}]

tmp2=Series[ellc2,{z,0,pos-1}]

tmp3=Series[ellc3,{z,0,pos-1}]

tmp4=Series[ell1,{z,0,pos-1}]

tmp5=Series[ell2,{z,0,pos-1}]

tmp6=Series[ell3,{z,0,pos-1}]

tmp7=Series[ell4,{z,0,pos-1}]

tmp8=Series[ell5,{z,0,pos-1}]

Series[D[S[z]^2ellc1,z],{z,0,2}]

Series[D[((z^2)/(S[z]))ellc2,z],{z,0,3}]

Print["Nice"];

Clear[p,h,a,T,S]

initp[z_]:=(z^2 p2+z^3 a0 h2+z^4 (((a1 h2)/2)-(2/3) (a0^2) p2)-z^5 ((1/3) (a0^3) h2+(2/3) a0 a1 p2+(1/3) p2 T3));

inith[z_]:=(z^2 h2-z^3 a0 p2-z^4 (((a1 p2)/2)+(2/3) (a0^2) h2)+z^5 ((1/3) (a0^3) p2-(2/3) a0 a1 h2-(1/3) h2 T3));

inita[z_]:=a0+z a1+(1/6)z^4 a0( (h2^2)+ (p2^2))+(1/6) z^5 a1 ( (h2^2)+ (p2^2));

initT[z_]:=1+z^3 T3+z^4 (((a1^2)/4)-((((h2^2)+(p2^2)))/3))+z^5 *0;

initS[z_]:=1-(1/6)z^4 ((h2^2)+(p2^2))+z^5 *0;

Check : Normal phase RN soln

p2=0;
h2=0;

eps=0.001;

zmax=2;

a0=1;
a1=-1 ;
T3=-1 ;

nsol=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[eps]==initp[eps],p'[eps]==initp'[eps],h[eps]==inith[eps],h'[eps]==inith'[eps],a[eps]==inita[eps],a'[eps]==inita'[eps],T[eps]==initT[eps],T'[eps]==initT'[eps],S[eps]==initS[eps],S'[eps]==initS'[eps]},{p,h,a,T,S},{z,eps,zmax}][[1]];

Plot[{T[z]/.nsol,1-z^3+(1/4)a1^2 z^4},{z,eps,zmax}]

Plot[S[z]/.nsol,{z,eps,zmax}]

Plot[p[z]/.nsol,{z,eps,zmax}]

Plot[h[z]/.nsol,{z,eps,zmax}]

Plot[{a[z]/.nsol,1+a1 z},{z,eps,zmax}]

Scalar soln

p2=+0.0001;
h2=0;

eps=0.001;

zmax=1.1;

a0=1;
a1=-1 ;
T3=-1 ;

nsol=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[eps]==initp[eps],p'[eps]==initp'[eps],h[eps]==inith[eps],h'[eps]==inith'[eps],a[eps]==inita[eps],a'[eps]==inita'[eps],T[eps]==initT[eps],T'[eps]==initT'[eps],S[eps]==initS[eps],S'[eps]==initS'[eps]},{p,h,a,T,S},{z,eps,zmax}][[1]];

Plot[{T[z]/.nsol,1-z^3+(1/4)a1^2 z^4},{z,eps,zmax}]

Plot[S[z]/.nsol,{z,eps,zmax}]

Plot[p[z]/.nsol,{z,eps,zmax}]

Plot[h[z]/.nsol,{z,eps,zmax}]

Plot[{a[z]/.nsol,1+a1 z},{z,eps,zmax}]

Print["\n\nHorizon expansion\n\n"]

Clear[a0,a1,T3,p2,h2]

pos=0;

T[z_]:=t1(1-z)+t2 (1-z)^2+t3(1-z)^3

S[z_]:=s0+s1(1-z)+s2 (1-z)^2+s3(1-z)^3

a[z_]:=a0+a1(1-z)+a2 (1-z)^2+a3(1-z)^3

p[z_]:=p0+p1(1-z)+p2 (1-z)^2+p3(1-z)^3

h[z_]:=h0+h1(1-z)+h2 (1-z)^2+h3(1-z)^3

tmp1=Simplify[Series[ell1,{z,1,pos}]]

tmp2=Simplify[Series[ell2,{z,1,pos}]]

tmp3=Simplify[Series[ell3,{z,1,pos-1}]]

tmp4=Simplify[Series[ell4,{z,1,pos-1}]]

tmp5=Simplify[Series[ell5,{z,1,pos-1}]]

tmp6=Simplify[Series[ellc1,{z,1,pos}]]

tmp7=Simplify[Series[ellc2,{z,1,pos}]]

tmp8=Simplify[Series[ellc3,{z,1,pos}]]

pos=0;

T[z_]:=t1(1-z)+t2 (1-z)^2+t3(1-z)^3

S[z_]:=s0+s1(1-z)+s2 (1-z)^2+s3(1-z)^3

a[z_]:=a1(1-z)+a2 (1-z)^2+a3(1-z)^3

p[z_]:=p0+p1(1-z)+p2 (1-z)^2+p3(1-z)^3

h[z_]:=h0+h1(1-z)+h2 (1-z)^2+h3(1-z)^3

tmp1=Simplify[Series[ell1,{z,1,pos}]]

tmp2=Simplify[Series[ell2,{z,1,pos}]]

tmp3=Simplify[Series[ell3,{z,1,pos-1}]]

tmp4=Simplify[Series[ell4,{z,1,pos-1}]]

tmp5=Simplify[Series[ell5,{z,1,pos-1}]]

tmp6=Simplify[Series[ellc1,{z,1,pos}]]

tmp7=Simplify[Series[ellc2,{z,1,pos}]]

tmp8=Simplify[Series[ellc3,{z,1,pos}]]

Simplify[Solve[{tmp1,tmp2,tmp7}==0,{h1,p1,s1}]]

pos=1;

T[z_]:=t1(1-z)+t2 (1-z)^2+t3(1-z)^3

S[z_]:=s0-((s0 ((a1^2)-4 (3+(h0^2)+(p0^2)-t1)))/(4 t1))(1-z)+s2 (1-z)^2+s3(1-z)^3

a[z_]:=a1(1-z)+a2 (1-z)^2+a3(1-z)^3

p[z_]:=p0+(-((a1 h0+2 p0)/t1))(1-z)+p2 (1-z)^2+p3(1-z)^3

h[z_]:=h0+(((-2 h0+a1 p0)/t1))(1-z)+h2 (1-z)^2+h3(1-z)^3

tmp1=Simplify[Series[ell1,{z,1,pos}]]

tmp2=Simplify[Series[ell2,{z,1,pos}]]

tmp3=Simplify[Series[ell3,{z,1,pos-1}]]

tmp4=Simplify[Series[ell4,{z,1,pos-1}]]

tmp5=Simplify[Series[ell5,{z,1,pos-1}]]

tmp6=Simplify[Series[ellc1,{z,1,pos}]]

tmp7=Simplify[Series[ellc2,{z,1,pos}]]

tmp8=Simplify[Series[ellc3,{z,1,pos}]]

Simplify[Solve[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp7}==0,{t2,p2,h2,s2,a2}]]

pos=1;

T[z_]:=t1(1-z)+((1/2) ((a1^2)-2 t1)) (1-z)^2+t3(1-z)^3

S[z_]:=s0-((s0 ((a1^2)-4 (3+(h0^2)+(p0^2)-t1)))/(4 t1))(1-z)+(-(((4+(a1^2)) ((h0^2)+(p0^2)) s0)/(4 (t1^2)))) (1-z)^2+s3(1-z)^3

a[z_]:=a1(1-z)+(((a1 ((a1^2)+4 (-3+t1)))/(4 t1))) (1-z)^2+a3(1-z)^3

p[z_]:=p0+(-((a1 h0+2 p0)/t1))(1-z)+((((a1^3) h0-4 (a1^2) p0+4 a1 h0 (7-2 t1)+8 p0 (4+(h0^2)+(p0^2)-2 t1))/(8 (t1^2)))) (1-z)^2+p3(1-z)^3

h[z_]:=h0+(((-2 h0+a1 p0)/t1))(1-z)+(-((-8 (h0^3)+a1 p0 (28+(a1^2)-8 t1)+4 h0 (-8+(a1^2)-2 (p0^2)+4 t1))/(8 (t1^2)))) (1-z)^2+h3(1-z)^3

tmp1=Simplify[Series[ell1,{z,1,pos}]]

tmp2=Simplify[Series[ell2,{z,1,pos}]]

tmp3=Simplify[Series[ell3,{z,1,pos-1}]]

tmp4=Simplify[Series[ell4,{z,1,pos-1}]]

tmp5=Simplify[Series[ell5,{z,1,pos-1}]]

tmp6=Simplify[Series[ellc1,{z,1,pos}]]

tmp7=Simplify[Series[ellc2,{z,1,pos}]]

tmp8=Simplify[Series[ellc3,{z,1,pos}]]

Clear[T,S,a,p,h]

initT[z_]:=t1(1-z)+((1/2) ((a1^2)-2 t1)) (1-z)^2

initS[z_]:=s0-((s0 ((a1^2)-4 (3+(h0^2)+(p0^2)-t1)))/(4 t1))(1-z)+(-(((4+(a1^2)) ((h0^2)+(p0^2)) s0)/(4 (t1^2)))) (1-z)^2

inita[z_]:=a1(1-z)+(((a1 ((a1^2)+4 (-3+t1)))/(4 t1))) (1-z)^2

initp[z_]:=p0+(-((a1 h0+2 p0)/t1))(1-z)+((((a1^3) h0-4 (a1^2) p0+4 a1 h0 (7-2 t1)+8 p0 (4+(h0^2)+(p0^2)-2 t1))/(8 (t1^2)))) (1-z)^2

inith[z_]:=h0+(((-2 h0+a1 p0)/t1))(1-z)+(-((-8 (h0^3)+a1 p0 (28+(a1^2)-8 t1)+4 h0 (-8+(a1^2)-2 (p0^2)+4 t1))/(8 (t1^2)))) (1-z)^2

Clear[p0,h0,a1,t1,s0]

p0=0;
h0=0;
a1=0;
t1=3;
s0=1;

zhoriz=1.00001;

eps=0.0001 ;

nsol=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[zhoriz]==initp[zhoriz],p'[zhoriz]==initp'[zhoriz],h[zhoriz]==inith[zhoriz],h'[zhoriz]==inith'[zhoriz],a[zhoriz]==inita[zhoriz],a'[zhoriz]==inita'[zhoriz],T[zhoriz]==initT[zhoriz],T'[zhoriz]==initT'[zhoriz],S[zhoriz]==initS[zhoriz],S'[zhoriz]==initS'[zhoriz]},{p,h,a,T,S},{z,eps,zhoriz}][[1]];

Plot[{T[z]/.nsol,1-z^3},{z,eps,zhoriz}]

Plot[S[z]/.nsol,{z,eps,zhoriz}]

Plot[p[z]/.nsol,{z,eps,zhoriz}]

Plot[h[z]/.nsol,{z,eps,zhoriz}]

Plot[{a[z]/.nsol},{z,eps,zhoriz}]

Clear[p0,a1,t1,s0]

p0=0.1;
h0=0.;
a1=1;
t1=2.8;
s0=1;

zhoriz=1.0000001;

eps=0.0001 ;

nsol=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[zhoriz]==initp[zhoriz],p'[zhoriz]==initp'[zhoriz],h[zhoriz]==inith[zhoriz],h'[zhoriz]==inith'[zhoriz],a[zhoriz]==inita[zhoriz],a'[zhoriz]==inita'[zhoriz],T[zhoriz]==initT[zhoriz],T'[zhoriz]==initT'[zhoriz],S[zhoriz]==initS[zhoriz],S'[zhoriz]==initS'[zhoriz]},{p,h,a,T,S},{z,eps,zhoriz}][[1]];

Plot[p[z]/.nsol,{z,eps,zhoriz}]

Plot[h[z]/.nsol,{z,eps,zhoriz}]

Clear[p0,h0,a1,t1,s0,inteq,zhoriz,eps];

inteq[p0in_Real,a1in_Real,t1in_Real,output_]:=(

p0=p0in;

h0=0;

a1=a1in;

t1=t1in;

s0=1;

eps=0.000001 ;

zhoriz=1-eps;

nsol=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[zhoriz]==initp[zhoriz],p'[zhoriz]==initp'[zhoriz],h[zhoriz]==inith[zhoriz],h'[zhoriz]==inith'[zhoriz],a[zhoriz]==inita[zhoriz],a'[zhoriz]==inita'[zhoriz],T[zhoriz]==initT[zhoriz],T'[zhoriz]==initT'[zhoriz],S[zhoriz]==initS[zhoriz],S'[zhoriz]==initS'[zhoriz]},{p,h,a,T,S},{z,eps,zhoriz}][[1]];

Clear[p0,h0,a1,t1,s0];

If[output==1,

Return[{h'[eps],T'[eps]}/.nsol];

];

If[output==0,

htab=Table[{z,h[z]/.nsol},{z,eps,0.001,0.0001}];

Ttab=Table[{z,T[z]-1/.nsol},{z,eps,0.001,0.0001}];

hfit=Fit[htab,{z,z^2},z];

Tfit=Fit[Ttab,z,z];

Return[{Coefficient[hfit,z],Coefficient[Tfit,z]}];

];

If[output==2,

ptab=Table[{z,p[z]/.nsol},{z,eps,0.001,0.0001}];

htab=Table[{z,h[z]/.nsol},{z,eps,0.001,0.0001}];

atab=Table[{z,a[z]/.nsol},{z,eps,0.001,0.0001}];

Ttab=Table[{z,T[z]-1/.nsol},{z,eps,0.001,0.0001}];

Stab=Table[{z,S[z]/.nsol},{z,eps,0.001,0.0001}];

pfit=Fit[ptab,{z,z^2},z];

hfit=Fit[htab,{z,z^2},z];

afit=Fit[atab,{1,z},z];

Tfit=Fit[Ttab,z,z];

Sfit=Fit[Stab,{1,z},z];

pfit2=Fit[ptab,{z^2},z];

hfit2=Fit[htab,{z^2},z];

Tfit2=Fit[Ttab,{z^3},z];

Return[{{pf==pfit,hf==hfit,af==afit,Tf==Tfit,Sf==Sfit},{pfit2,hfit2,Tfit2}}];

];

);

FindRoot[inteq[x,3.,y,1]=={0,0},{x,2.2},{y,3.}]

solxy=FindRoot[inteq[x,3.,y,1]=={0,0},{x,2.2},{y,3}]

solxy2=FindRoot[inteq[x,3.,y,0]=={0,0},{x,(x/.solxy)},{y,(y/.solxy)}]

inteq[x,3.,y,1]/.solxy2

inteq[x,3.,y,2]/.solxy2

Plot[T[z]/.nsol,{z,eps,zhoriz},PlotRange->All]

Plot[a[z]/.nsol,{z,eps,zhoriz},PlotRange->All]

Plot[{p[z],h[z]}/z/.nsol,{z,eps,zhoriz},PlotRange->All]

sol=FullSimplify[Solve[{ell1,ell2,ell3,ell4,ell5}==0,{p''[z],h''[z],a''[z],T''[z],S''[z]}][[1]]]

sol2=sol/.{p->fp,h->fh,a->fa,T->fT,S->fS};

p[z_]:=fp[z]Cos[α]-fh[z] Sin[α];

h[z_]:=fh[z]Cos[α]+fp[z]Sin[α];

a[z_]:=fa[z];

T[z_]:=fT[z];

S[z_]:=kk fS[z];

Simplify[ell1/.sol2]

Simplify[ell2/.sol2]

Simplify[ell3/.sol2]

Simplify[ell4/.sol2]

Simplify[ell5/.sol2]

Numerics

Clear[p,h,a,T,S]

ell1=a[z] (h'[z]+h[z] (-(1/z)+((S'[z])/(S[z]))))+(1/2) (-((2 p[z])/(z^2))+h[z] a'[z]+((2 T[z] p'[z])/z)-((2 T[z] p'[z] S'[z])/(S[z]))-p'[z] T'[z]-T[z] p''[z]);

ell2=-(1/(2 (z^2) S[z]))(2 (z^2) T[z] h'[z] S'[z]+2 z a[z] (z S[z] p'[z]+p[z] (-S[z]+z S'[z]))+S[z] (2 h[z]+z (z p[z] a'[z]-2 T[z] h'[z]+z h'[z] T'[z]+z T[z] h''[z])));

ellc1=((2 (a[z] (((h[z])^2)+((p[z])^2))+T[z] (p[z] h'[z]-h[z] p'[z])))/(z^2));

ellc2=-((4 (z^2) T[z] ((S'[z])^2)+4 z S[z] S'[z] (-4 T[z]+z T'[z])-((S[z])^2) (12+4 ((h[z])^2)+4 ((p[z])^2)-(z^4) ((a'[z])^2)+2 T[z] (-6+(z^2) ((h'[z])^2)+(z^2) ((p'[z])^2))+4 z T'[z]))/(8 (z^2) S[z]));

ell3=-((2 z T[z] ((S'[z])^2)+((S[z])^2) ((z^3) ((a'[z])^2)+2 (z a[z] (-p[z] h'[z]+h[z] p'[z])-z T[z] (((h'[z])^2)+((p'[z])^2))+T'[z]))-2 S[z] T[z] (2 S'[z]+z S''[z]))/(z ((S[z])^2)))+T''[z];

ell4=(1/2) S[z] (((h'[z])^2)+((p'[z])^2))+S''[z];

ell5=((2 p[z] h'[z])/(z^2))-((2 h[z] p'[z])/(z^2))+((2 a'[z] S'[z])/(S[z]))+a''[z];

ellc3=(1/(8 (z^2) S[z]))(4 (z^2) ((a[z])^2) (((h[z])^2)+((p[z])^2)) ((S[z])^2)+4 (z^2) a[z] ((S[z])^2) T[z] (p[z] h'[z]-h[z] p'[z])+T[z] (-4 (z^2) T[z] ((S'[z])^2)-4 z S[z] S'[z] (-4 T[z]+z T'[z])+((S[z])^2) (12+4 ((h[z])^2)+4 ((p[z])^2)-12 T[z]-(z^4) ((a'[z])^2)+2 (z^2) T[z] ((h'[z])^2)+2 (z^2) T[z] ((p'[z])^2)+4 z T'[z])));

tmp=Simplify[Solve[{ellc1==0,ellc2==0},{T'[z],h'[z]}]]

Expand[tmp]

Simplify[ellc3/.tmp]

tmp=Simplify[Solve[{ell1==0,ell4==0,ell5==0,ellc1==0,ellc2==0},{p''[z],a''[z],S''[z],T'[z],h'[z]}]]

eqlist=Thread[{p''[z],a''[z],S''[z],T'[z],h'[z]}==({p''[z],a''[z],S''[z],T'[z],h'[z]}/.tmp[[1]])]

eqlist={p''[z]==-(2 (z^2) ((a[z])^2) (((h[z])^2)+((p[z])^2)) S[z] (z ((h[z])^2) S[z] p'[z]+z ((p[z])^2) S[z] p'[z]-4 p[z] (S[z]-z S'[z]))-4 z a[z] h[z] T[z] ((z^2) ((h[z])^2) ((S[z])^2) ((p'[z])^2)+2 z p[z] S[z] p'[z] (-S[z]+z S'[z])+((p[z])^2) (((S[z])^2) (2+(z^2) ((p'[z])^2))-4 z S[z] S'[z]+2 (z^2) ((S'[z])^2)))+T[z] (4 z ((p[z])^4) ((S[z])^2) p'[z]+2 (z^3) ((h[z])^2) ((S[z])^2) T[z] ((p'[z])^3)-8 ((p[z])^3) S[z] (S[z]-z S'[z])+z ((p[z])^2) (4 ((h[z])^2) ((S[z])^2) p'[z]+4 z h[z] S[z] a'[z] (S[z]-z S'[z])+p'[z] (((S[z])^2) (12-(z^4) ((a'[z])^2)+2 T[z] (-2+(z^2) ((p'[z])^2)))+4 (z^2) T[z] ((S'[z])^2)))))/(4 (z^2) ((p[z])^2) S[z] ((T[z])^2) (-S[z]+z S'[z])),a''[z]==((2 (a[z] (((h[z])^2)+((p[z])^2)) S[z]-(z^2) T[z] a'[z] S'[z]))/((z^2) S[z] T[z])),S''[z]==-(((((h[z])^2)+((p[z])^2)) S[z] (((a[z])^2) (((h[z])^2)+((p[z])^2))-2 a[z] h[z] T[z] p'[z]+((T[z])^2) ((p'[z])^2)))/(2 ((p[z])^2) ((T[z])^2))),T'[z]==(2 (z^2) ((a[z])^2) (((((h[z])^2)+((p[z])^2)))^2) ((S[z])^2)-4 (z^2) a[z] h[z] (((h[z])^2)+((p[z])^2)) ((S[z])^2) T[z] p'[z]+T[z] (4 ((p[z])^4) ((S[z])^2)+2 (z^2) ((h[z])^2) ((S[z])^2) T[z] ((p'[z])^2)+((p[z])^2) (((S[z])^2) (12+4 ((h[z])^2)-(z^4) ((a'[z])^2)+2 T[z] (-6+(z^2) ((p'[z])^2)))+16 z S[z] T[z] S'[z]-4 (z^2) T[z] ((S'[z])^2))))/(4 z ((p[z])^2) S[z] T[z] (-S[z]+z S'[z])),h'[z]==((-a[z] (((h[z])^2)+((p[z])^2))+h[z] T[z] p'[z])/(p[z] T[z]))};

repl1={p2->asymp2,h2->asymh2,a0->asyma0,a1->asyma1,T3->asymT3};

repl2={t1->horizt1,a1->horiza1,p0->horizp0,h0->horizh0,s0->horizs0};

asymp[z_]:=(z^2 p2+z^3 a0 h2+z^4 (((a1 h2)/2)-(2/3) (a0^2) p2)-z^5 ((1/3) (a0^3) h2+(2/3) a0 a1 p2+(1/3) p2 T3))/.repl1;

asymh[z_]:=(z^2 h2-z^3 a0 p2-z^4 (((a1 p2)/2)+(2/3) (a0^2) h2)+z^5 ((1/3) (a0^3) p2-(2/3) a0 a1 h2-(1/3) h2 T3))/.repl1;

asyma[z_]:=(a0+z a1+(1/6)z^4 a0( (h2^2)+ (p2^2))+(1/6) z^5 a1 ( (h2^2)+ (p2^2)))/.repl1;

asymT[z_]:=(1+z^3 T3+z^4 (((a1^2)/4)-((((h2^2)+(p2^2)))/3)))/.repl1;

asymS[z_]:=(1-(1/6)z^4 ((h2^2)+(p2^2)))/.repl1;

horizT[z_]:=(t1(1-z)+((1/2) ((a1^2)-2 t1)) (1-z)^2)/.repl2

horizS[z_]:=(s0-((s0 ((a1^2)-4 (3+(h0^2)+(p0^2)-t1)))/(4 t1))(1-z)+(-(((4+(a1^2)) ((h0^2)+(p0^2)) s0)/(4 (t1^2)))) (1-z)^2)/.repl2

horiza[z_]:=(a1(1-z)+(((a1 ((a1^2)+4 (-3+t1)))/(4 t1))) (1-z)^2)/.repl2

horizp[z_]:=(p0+(-((a1 h0+2 p0)/t1))(1-z)+((((a1^3) h0-4 (a1^2) p0+4 a1 h0 (7-2 t1)+8 p0 (4+(h0^2)+(p0^2)-2 t1))/(8 (t1^2)))) (1-z)^2)/.repl2

horizh[z_]:=(h0+(((-2 h0+a1 p0)/t1))(1-z)+(-((-8 (h0^3)+a1 p0 (28+(a1^2)-8 t1)+4 h0 (-8+(a1^2)-2 (p0^2)+4 t1))/(8 (t1^2)))) (1-z)^2)/.repl2

Clear[inteq,zhoriz,eps];

inteq[horizp0in_Real,horiza1in_Real,horizt1in_Real,asymp2in_Real,asymh2in_Real,asyma0in_Real,asyma1in_Real,asymT3in_Real,output_]:=(

horizp0=horizp0in;

horizh0=0;

horiza1=horiza1in;

horizt1=horizt1in;

horizs0=1;

asymp2=asymp2in;

asymh2=asymh2in;

asyma0=asyma0in;

asyma1=asyma1in;

asymT3=asymT3in;


eps=0.01 ;

zmid = 0.5 ;

zhoriz=1-eps;

(*nsol1=NDSolve[{ell1==0,ell2==0,ell4==0,ellc1==0,ellc2==0,p[eps]==asymp[eps],p'[eps]==asymp'[eps],h[eps]==asymh[eps],h'[eps]==asymh'[eps],T[eps]==asymT[eps],S[eps]==asymS[eps],S'[eps]==asymS'[eps]},{p,h,a,T,S},{z,eps,zhoriz}][[1]];
*)


nsol2=NDSolve[Join[eqlist,{p[zhoriz]==horizp[zhoriz],p'[zhoriz]==horizp'[zhoriz],a[zhoriz]==horiza[zhoriz],a'[zhoriz]==horiza'[zhoriz],S[zhoriz]==horizS[zhoriz],S'[zhoriz]==horizS'[zhoriz],T[zhoriz]==horizT[zhoriz],h[zhoriz]==horizh[zhoriz]}],{p,h,a,T,S},{z,eps,zhoriz}][[1]];

Clear[horizp0,horizh0,horiza1,horizt1,horizs0,asymp2,asymh2,asyma0,asyma1,asymT3];

If[output==1,

Return[{h'[eps],T'[eps]}/.nsol2];

];


If[output==0,

ptab=Table[{z,p[z]/.nsol2},{z,eps,0.001,0.0001}];

htab=Table[{z,h[z]/.nsol2},{z,eps,0.001,0.0001}];

atab=Table[{z,a[z]/.nsol2},{z,eps,0.001,0.0001}];

Ttab=Table[{z,T[z]-1/.nsol2},{z,eps,0.001,0.0001}];

Stab=Table[{z,S[z]/.nsol2},{z,eps,0.001,0.0001}];

pfit=Fit[ptab,{z,z^2},z];

hfit=Fit[htab,{z,z^2},z];

afit=Fit[atab,{1,z},z];

Tfit=Fit[Ttab,z,z];

Sfit=Fit[Stab,{1,z},z];

pfit2=Fit[ptab,{z^2},z];

hfit2=Fit[htab,{z^2},z];

Tfit2=Fit[Ttab,{z^3},z];

Return[{{pf==pfit,hf==hfit,af==afit,Tf==Tfit,Sf==Sfit},{pfit2,hfit2,Tfit2}}];

];

);

inteq[2.2,3.,0.16,0.,0.,0.,0.,0.,1]

Clear[inteq,zhoriz,eps];

inteq[horizp0in_Real,horiza1in_Real,horizt1in_Real,asymp2in_Real,asymh2in_Real,asyma0in_Real,asyma1in_Real,asymT3in_Real,output_]:=(

horizp0=horizp0in;

horizh0=0;

horiza1=horiza1in;

horizt1=horizt1in;

horizs0=1;

asymp2=asymp2in;

asymh2=asymh2in;

asyma0=asyma0in;

asyma1=asyma1in;

asymT3=asymT3in;


eps=0.000001 ;

eps2=0.000001 ;

zmid = 0.2;

zhoriz=1-eps2;

nsol1=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[eps]==asymp[eps],p'[eps]==asymp'[eps],h[eps]==asymh[eps],h'[eps]==asymh'[eps],a[eps]==asyma[eps],a'[eps]==asyma'[eps],T[eps]==asymT[eps],T'[eps]==asymT'[eps],S[eps]==asymS[eps],S'[eps]==asymS'[eps]},{p,h,a,T,S},{z,eps,zhoriz}][[1]];

zcut1=nsol1[[1,2,1,1,2]];

If[zcut1<zmid,Print["Error"]];

nsol2=NDSolve[{ell1==0,ell2==0,ell3==0,ell4==0,ell5==0,p[zhoriz]==horizp[zhoriz],p'[zhoriz]==horizp'[zhoriz],h[zhoriz]==horizh[zhoriz],h'[zhoriz]==horizh'[zhoriz],a[zhoriz]==horiza[zhoriz],a'[zhoriz]==horiza'[zhoriz],T[zhoriz]==horizT[zhoriz],T'[zhoriz]==horizT'[zhoriz],S[zhoriz]==horizS[zhoriz],S'[zhoriz]==horizS'[zhoriz]},{p,h,a,T,S},{z,eps,zhoriz}][[1]];

zcut2=nsol2[[1,2,1,1,1]];

If[zcut2>zmid,Print["Error"]];

Clear[horizp0,horizh0,horiza1,horizt1,horizs0,asymp2,asymh2,asyma0,asyma1,asymT3];

If[output==1,

(*Return[{h'[eps],T'[eps]}/.nsol2];*)

htab=Table[{z,h[z]/.nsol2},{z,20eps,100eps,10eps}];

Ttab=Table[{z,T[z]-1/.nsol2},{z,20eps,100eps,10eps}];

hfit=Fit[htab,{z,z^2},z];

Tfit=Fit[Ttab,{z,z^2},z];

Return[{Coefficient[hfit,z],Coefficient[Tfit,z]}];


];

If[output==2,

list={(p[zmid]/.nsol1)-(p[zmid]/.nsol2),(h[zmid]/.nsol1)-(h[zmid]/.nsol2),(a[zmid]/.nsol1)-(a[zmid]/.nsol2),(T[zmid]/.nsol1)-(T[zmid]/.nsol2)};

dlist={(p'[zmid]/.nsol1)-(p'[zmid]/.nsol2),(h'[zmid]/.nsol1)-(h'[zmid]/.nsol2),(a'[zmid]/.nsol1)-(a'[zmid]/.nsol2),(T'[zmid]/.nsol1)-(T'[zmid]/.nsol2),(((S'[zmid])/(S[zmid]))/.nsol1)-(((S'[zmid])/(S[zmid]))/.nsol2)};

Return[Join[list,{dlist[[1]],dlist[[2]],dlist[[3]]}]];

];


If[output==0,

ptab=Table[{z,p[z]/.nsol2},{z,20eps,100eps,10eps}];

htab=Table[{z,h[z]/.nsol2},{z,20eps,100eps,10eps}];

atab=Table[{z,a[z]/.nsol2},{z,20eps,100eps,10eps}];

Ttab=Table[{z,T[z]-1/.nsol2},{z,20eps,100eps,10eps}];

Stab=Table[{z,S[z]/.nsol2},{z,20eps,100eps,10eps}];

pfit=Fit[ptab,{z,z^2},z];

hfit=Fit[htab,{z,z^2},z];

afit=Fit[atab,{1,z},z];

Tfit=Fit[Ttab,{z,z^2,z^3},z];

Sfit=Fit[Stab,{1,z},z];

pfit2=Fit[ptab,{z^2},z];

hfit2=Fit[htab,{z^2},z];

Tfit2=Fit[Ttab,{z^3},z];

Return[{{pf==pfit,hf==hfit,af==afit,Tf==Tfit,Sf==Sfit},{pfit2,hfit2,Tfit2}}];

];

);

inteq[2.2,3.,0.16,0.,0.,0.,0.,0.,1]

solxy=FindRoot[inteq[x,3.,y,0.,0.,0.,0.,0.,1]=={0,0},{x,2.2},{y,0.2}]

inteq[x,3.,y,0.,0.,0.,0.,0.,0]/.solxy

inteq[x,3.,y,-1.382662769389146`,-0.26907561405383607` ,3.2494537744400596`,-3.454184133737518`  ,-3.73,0]/.solxy

inteq[x,3.,y,-1.382662769389146`,-0.26907561405383607` ,3.2494537744400596`,-3.454184133737518`  ,-3.73,2]/.solxy

solall=FindRoot[inteq[x1,3.,x3,y1,y2,y3,y4,y5,2]==0,{x1,2.1},{x3,0.13},{y1,-1.4},{y2,-0.27},{y3,3.},{y4,-3.5},{y5,-3.73}]

inteq[x1,3.,x3,y1,y2,y3,y4,y5,2]/.solall

plt1=Plot[p[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[p[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[h[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[h[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[a[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[a[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[T[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[T[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[((S[z])/(S[zmid]))/.nsol1,{z,eps,zmid}];
plt2=Plot[((S[z])/(S[zmid]))/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

eqlist

Clear[inteq,zhoriz,eps];

inteq[horizp0in_Real,horiza1in_Real,horizt1in_Real,horizs0in_Real,asymp2in_Real,asymh2in_Real,asyma0in_Real,asyma1in_Real,asymT3in_Real,output_]:=(

horizp0=horizp0in;

horizh0=0;

horiza1=horiza1in;

horizt1=horizt1in;

horizs0=horizs0in;

asymp2=asymp2in;

asymh2=asymh2in;

asyma0=asyma0in;

asyma1=asyma1in;

asymT3=asymT3in;


eps=0.0001 ;

eps2=0.000001 ;

zmid = 0.2;

zhoriz=1-eps2;

nsol1=NDSolve[Join[eqlist,{p[eps]==asymp[eps],p'[eps]==asymp'[eps],h[eps]==asymh[eps],a[eps]==asyma[eps],a'[eps]==asyma'[eps],T[eps]==asymT[eps],S[eps]==asymS[eps],S'[eps]==asymS'[eps]}],{p,h,a,T,S},{z,eps,zhoriz}][[1]];

zcut1=nsol1[[1,2,1,1,2]];

If[zcut1<zmid,Print["Error"]];

nsol2=NDSolve[Join[eqlist,{p[zhoriz]==horizp[zhoriz],p'[zhoriz]==horizp'[zhoriz],h[zhoriz]==horizh[zhoriz],a[zhoriz]==horiza[zhoriz],a'[zhoriz]==horiza'[zhoriz],T[zhoriz]==horizT[zhoriz],S[zhoriz]==horizS[zhoriz],S'[zhoriz]==horizS'[zhoriz]}],{p,h,a,T,S},{z,eps,zhoriz}][[1]];

zcut2=nsol2[[1,2,1,1,1]];

If[zcut2>zmid,Print["Error"]];

Clear[horizp0,horizh0,horiza1,horizt1,horizs0,asymp2,asymh2,asyma0,asyma1,asymT3];

If[output==1,

(*Return[{h'[eps],T'[eps]}/.nsol2];*)

htab=Table[{z,h[z]/.nsol2},{z,20eps,100eps,10eps}];

Ttab=Table[{z,T[z]-1/.nsol2},{z,20eps,100eps,10eps}];

hfit=Fit[htab,{z,z^2},z];

Tfit=Fit[Ttab,{z,z^2},z];

Return[{Coefficient[hfit,z],Coefficient[Tfit,z]}];


];

If[output==2,

list={(p[zmid]/.nsol1)-(p[zmid]/.nsol2),(h[zmid]/.nsol1)-(h[zmid]/.nsol2),(a[zmid]/.nsol1)-(a[zmid]/.nsol2),(T[zmid]/.nsol1)-(T[zmid]/.nsol2),(S[zmid]/.nsol1)-(S[zmid]/.nsol2)};

dlist={(p'[zmid]/.nsol1)-(p'[zmid]/.nsol2),(a'[zmid]/.nsol1)-(a'[zmid]/.nsol2),(((S'[zmid])/(S[zmid]))/.nsol1)-(((S'[zmid])/(S[zmid]))/.nsol2)};

Return[Join[list,{dlist[[1]],dlist[[2]],dlist[[3]]}]];

];


If[output==0,

ptab=Table[{z,p[z]/.nsol2},{z,20eps,100eps,10eps}];

htab=Table[{z,h[z]/.nsol2},{z,20eps,100eps,10eps}];

atab=Table[{z,a[z]/.nsol2},{z,20eps,100eps,10eps}];

Ttab=Table[{z,T[z]-1/.nsol2},{z,20eps,100eps,10eps}];

Stab=Table[{z,S[z]/.nsol2},{z,20eps,100eps,10eps}];

pfit=Fit[ptab,{z,z^2},z];

hfit=Fit[htab,{z,z^2},z];

afit=Fit[atab,{1,z},z];

Tfit=Fit[Ttab,{z,z^2,z^3},z];

Sfit=Fit[Stab,{1,z},z];

pfit2=Fit[ptab,{z^2},z];

hfit2=Fit[htab,{z^2},z];

Tfit2=Fit[Ttab,{z^3},z];

Return[{{pf==pfit,hf==hfit,af==afit,Tf==Tfit,Sf==Sfit},{pfit2,hfit2,Tfit2}}];

];

);

solxy=FindRoot[inteq[x,3.,y,1.,-1.382662769389146`,-0.26907561405383607` ,3.2494537744400596`,-3.454184133737518`  ,-3.73,1]=={0,0},{x,2.2},{y,0.2}]

inteq[x,3.,y,1.,-1.382662769389146`,-0.26907561405383607` ,3.2494537744400596`,-3.454184133737518`  ,-3.73,0]/.solxy

inteq[x,3.,y,0.3,-1.382662769389146`,-0.26907561405383607` ,3.2494537744400596`,-3.454184133737518`  ,-3.73,2]/.solxy

solall=FindRoot[inteq[x1,3.,x3,x4,y1,y2,y3,y4,y5,2]==0,{x1,2.27},{x3,0.12},{x4,1.},{y1,-1.38},{y2,-0.27},{y3,3.3},{y4,-3.4},{y5,-3.73}]

inteq[x1,3.,x3,y1,y2,y3,y4,y5,2]/.solall

plt1=Plot[p[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[p[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[h[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[h[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[a[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[a[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[T[z]/.nsol1,{z,eps,zmid}];
plt2=Plot[T[z]/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

plt1=Plot[((S[z])/(S[zmid]))/.nsol1,{z,eps,zmid}];
plt2=Plot[((S[z])/(S[zmid]))/.nsol2,{z,zmid,zhoriz}];
Show[plt1,plt2,PlotRange->All]

vec={2.27,0.12,0.3,-1.38,-0.27,3.3,-3.4,-3.73};

Do[

mat={};

vec0=vec;

eq0=inteq[vec0[[1]],3.,vec0[[2]],vec0[[3]],vec0[[4]],vec0[[5]],vec0[[6]],vec0[[7]],vec0[[8]],2];

DX = 10^-3;

Do[

vec0[[jj]]-=DX ;

eq1=inteq[vec0[[1]],3.,vec0[[2]],vec0[[3]],vec0[[4]],vec0[[5]],vec0[[6]],vec0[[7]],vec0[[8]],2];

vec0[[jj]]+=2*DX ;

eq2=inteq[vec0[[1]],3.,vec0[[2]],vec0[[3]],vec0[[4]],vec0[[5]],vec0[[6]],vec0[[7]],vec0[[8]],2];

AppendTo[mat,(eq2-eq1)/(2*DX)];

,{jj,1,8}];

dvec=Inverse[mat].eq0;

Print["ii = ",ii," -> dvec = ",dvec];

vec = vec-0.01*dvec;

,{ii,1,1}];

"inteq[horizp0in_Real,horiza1in_Real,horizt1in_Real,horizs0in_Real,asymp2in_Real,asymh2in_Real,asyma0in_Real,asyma1in_Real,asymT3in_Real,output_]"

vec={2.27,0.12,0.3,-1.38,-0.27,3.3,-3.4,-3.73};

Series[1-z^3+((a1^2)/4) z^4,{z,1,1}]

Series[1-z,{z,1,1}]

vec={0,0.12,0.3,-1.38,-0.27,3.3,-3.4,-3.73};

mat={};

eq0=inteq[vec[[1]],3.,vec[[2]],vec[[3]],vec[[4]],vec[[5]],vec[[6]],vec[[7]],vec[[8]],2];

DX = 10^-4;

Do[

vec0=vec;

vec0[[jj]]-=DX ;

eq1=inteq[vec0[[1]],3.,vec0[[2]],vec0[[3]],vec0[[4]],vec0[[5]],vec0[[6]],vec0[[7]],vec0[[8]],2];

vec0[[jj]]+=2*DX ;

eq2=inteq[vec0[[1]],3.,vec0[[2]],vec0[[3]],vec0[[4]],vec0[[5]],vec0[[6]],vec0[[7]],vec0[[8]],2];

AppendTo[mat,(eq2-eq1)/(2*DX)];

,{jj,1,8}];

dvec=Inverse[mat].eq0;

vec = vec-0.001*dvec;

Print["vec = ",vec,"  ->  eq = ",eq0];

vec

mat

inteq[x1,3.,x3,y1,y2,y3,y4,y5,2]==0,{x1,2.27},{x3,0.12},{y1,-1.38},{y2,-0.27},{y3,3.3},{y4,-3.4},{y5,-3.73}]