Previous notebook

SetDirectory[NotebookDirectory[]];

ClearRepl

Aseries[z_,t_,ρ_]=(r^2)+2 r ξ[t,ρ]+((ξ[t,ρ])^2)+(((Subscript[a, 4])[t,ρ])/(r^2))-2 Derivative[1,0][ξ][t,ρ]+((-2 ξ[t,ρ] (Subscript[a, 4])[t,ρ]+(1/2) Derivative[1,0][(Subscript[a, 4])][t,ρ])/(r^3))+(1/(60 (r^4) (t^6) ρ))(-8 ρ+t (-8 (t^2) ρ ((ξ[t,ρ])^3)+12 ((ξ[t,ρ])^2) (-2 t ρ+(t^3) (t (15 t ρ (Subscript[a, 4])[t,ρ]+Derivative[0,1][ξ][t,ρ]+ρ Derivative[0,2][ξ][t,ρ])+2 ρ Derivative[1,0][ξ][t,ρ]))+6 ξ[t,ρ] (2 (t^3) Derivative[0,1][ξ][t,ρ]+4 (t^4) ρ ((Derivative[0,1][ξ][t,ρ])^2)+ρ (-4+2 (t^3) Derivative[0,2][ξ][t,ρ]+4 (t^2) Derivative[1,0][ξ][t,ρ]-15 (t^5) Derivative[1,0][(Subscript[a, 4])][t,ρ]))+t (12 (t^3) (Subscript[f, 4])[t,ρ]+4 t Derivative[0,1][ξ][t,ρ]+12 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)-3 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-12 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]+6 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]+12 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+4 t ρ Derivative[0,2][ξ][t,ρ]-3 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-6 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]+8 ρ Derivative[1,0][ξ][t,ρ]+3 (t^3) ρ Derivative[1,0][(Subscript[a, 4])][t,ρ]-6 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]-6 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]+9 (t^4) ρ Derivative[2,0][(Subscript[a, 4])][t,ρ])))+(1/(270 (r^5) (t^7) (ρ^2)))(140 (ρ^2)+2208 (t^4) (ρ^2) ((ξ[t,ρ])^4)+1440 (t^5) (ρ^2) ((ξ[t,ρ])^5)+12 (t^4) (ρ^2) (Subscript[a, 4])[t,ρ]-8 (t^3) ρ ((ξ[t,ρ])^3) (-223 ρ+27 (t^2) (t (5 t ρ (Subscript[a, 4])[t,ρ]+Derivative[0,1][ξ][t,ρ]+ρ Derivative[0,2][ξ][t,ρ])+2 ρ Derivative[1,0][ξ][t,ρ]))-12 (t^2) ρ ((ξ[t,ρ])^2) (-80 ρ+90 (t^4) ρ (Subscript[a, 4])[t,ρ]+90 (t^4) ρ (Subscript[b, 4])[t,ρ]+3 (t^2) (30 (t^2) ρ (Subscript[c, 4])[t,ρ]+30 (t^3) (Subscript[f, 4])[t,ρ]+9 t Derivative[0,1][ξ][t,ρ]+12 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+30 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+7 t ρ Derivative[0,2][ξ][t,ρ]+18 ρ Derivative[1,0][ξ][t,ρ]-(t^2) Derivative[1,1][ξ][t,ρ]-(t^2) ρ Derivative[1,2][ξ][t,ρ]-2 t ρ Derivative[2,0][ξ][t,ρ]))+2 t ρ ξ[t,ρ] (232 ρ+9 (t^2) (4 (t^2) ρ (Subscript[a, 4])[t,ρ]+4 (t^2) ρ (Subscript[b, 4])[t,ρ]+4 (t^2) ρ (Subscript[c, 4])[t,ρ]-8 (t^3) (Subscript[f, 4])[t,ρ]-12 t Derivative[0,1][ξ][t,ρ]-16 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+3 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]+12 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]-6 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]-8 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]-8 t ρ Derivative[0,2][ξ][t,ρ]+3 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]+6 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-28 ρ Derivative[1,0][ξ][t,ρ]+4 (t^3) Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+4 (t^3) ρ Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+8 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+6 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]+6 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]+2 (t^2) Derivative[1,1][ξ][t,ρ]+8 (t^3) ρ Derivative[0,1][ξ][t,ρ] Derivative[1,1][ξ][t,ρ]+2 (t^2) ρ Derivative[1,2][ξ][t,ρ]+4 t ρ Derivative[2,0][ξ][t,ρ]-9 (t^4) ρ Derivative[2,0][(Subscript[a, 4])][t,ρ]))-3 (t^2) (-4 (t^2) (3 (t^2)+(ρ^2)) (Subscript[b, 4])[t,ρ]+4 (3 (t^4)-(t^2) (ρ^2)) (Subscript[c, 4])[t,ρ]+ρ (20 (t^3) (Subscript[f, 4])[t,ρ]+3 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-12 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]+8 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+12 t ρ Derivative[0,2][ξ][t,ρ]+56 ρ Derivative[1,0][ξ][t,ρ]-12 (t^3) ρ Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]-24 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)-24 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2) (-1+(t^2) Derivative[1,0][ξ][t,ρ])-6 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]-6 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]-4 (t^2) Derivative[1,1][ξ][t,ρ]-4 t Derivative[0,1][ξ][t,ρ] (-5+3 (t^2) (Derivative[1,0][ξ][t,ρ]+2 ρ Derivative[1,1][ξ][t,ρ]))+3 (t^5) Derivative[1,1][(Subscript[a, 4])][t,ρ]+12 (t^5) Derivative[1,1][(Subscript[b, 4])][t,ρ]-6 (t^5) Derivative[1,1][(Subscript[c, 4])][t,ρ]-12 (t^4) ρ Derivative[1,1][(Subscript[f, 4])][t,ρ]-4 (t^2) ρ Derivative[1,2][ξ][t,ρ]+3 (t^5) ρ Derivative[1,2][(Subscript[a, 4])][t,ρ]+6 (t^5) ρ Derivative[1,2][(Subscript[b, 4])][t,ρ]-8 t ρ Derivative[2,0][ξ][t,ρ]-3 (t^4) ρ Derivative[2,0][(Subscript[a, 4])][t,ρ]+6 (t^4) ρ Derivative[2,0][(Subscript[b, 4])][t,ρ]+6 (t^4) ρ Derivative[2,0][(Subscript[c, 4])][t,ρ]-3 (t^5) ρ Derivative[3,0][(Subscript[a, 4])][t,ρ])))/.r->1/z;
Bseries[z_,t_,ρ_]=-(2/(3 r t))-(2/3) Log[t]+((1+2 t ξ[t,ρ])/(3 (r^2) (t^2)))-((2 (1+3 t ξ[t,ρ]+3 (t^2) ((ξ[t,ρ])^2)))/(9 (r^3) (t^3)))+(((Subscript[b, 4])[t,ρ])/(r^4))+(((8/(15 (t^5)))+((2 ξ[t,ρ])/(t^4))+((10 ((ξ[t,ρ])^2))/(3 (t^3)))+((10 ((ξ[t,ρ])^3))/(3 (t^2)))+((2 ((ξ[t,ρ])^4))/t)+((2 (Subscript[a, 4])[t,ρ])/(15 t))-4 ξ[t,ρ] (Subscript[b, 4])[t,ρ]+((2 (Subscript[f, 4])[t,ρ])/(15 ρ))-(4/15) Derivative[0,1][(Subscript[f, 4])][t,ρ]-((2 Derivative[1,0][ξ][t,ρ])/(3 (t^3)))-((2 ξ[t,ρ] Derivative[1,0][ξ][t,ρ])/(t^2))-((2 ((ξ[t,ρ])^2) Derivative[1,0][ξ][t,ρ])/t)+Derivative[1,0][(Subscript[b, 4])][t,ρ])/(r^5))+(1/(540 (r^6) (t^6) ρ))(-960 ρ+t (-6300 (t^3) ρ ((ξ[t,ρ])^4)-3240 (t^4) ρ ((ξ[t,ρ])^5)+90 (t^2) ρ ((ξ[t,ρ])^3) (-91+60 (t^2) Derivative[1,0][ξ][t,ρ])+30 ((ξ[t,ρ])^2) (180 (t^5) ρ (Subscript[b, 4])[t,ρ]+5 (t^4) Derivative[0,1][ξ][t,ρ]-t ρ (240+(t^2) (t Derivative[0,2][ξ][t,ρ]-219 Derivative[1,0][ξ][t,ρ]+21 t Derivative[2,0][ξ][t,ρ])))+30 ξ[t,ρ] (-129 ρ-12 (t^4) ρ (Subscript[a, 4])[t,ρ]-(t^2) (12 (t^3) (Subscript[f, 4])[t,ρ]-5 t Derivative[0,1][ξ][t,ρ]+2 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+ρ (-141 Derivative[1,0][ξ][t,ρ]+t (Derivative[0,2][ξ][t,ρ]+6 t (-4 t Derivative[0,1][(Subscript[f, 4])][t,ρ]+7 ((Derivative[1,0][ξ][t,ρ])^2)+15 t Derivative[1,0][(Subscript[b, 4])][t,ρ])+21 Derivative[2,0][ξ][t,ρ]))))+t (-72 (t^2) ρ (Subscript[a, 4])[t,ρ]+60 (t^3) (Subscript[f, 4])[t,ρ]+50 t Derivative[0,1][ξ][t,ρ]-30 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+12 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-15 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]-60 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]-30 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]-10 t ρ Derivative[0,2][ξ][t,ρ]-24 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]+15 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]+1230 ρ Derivative[1,0][ξ][t,ρ]-630 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+72 (t^3) ρ Derivative[1,0][(Subscript[a, 4])][t,ρ]+45 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]+90 (t^4) Derivative[1,0][(Subscript[f, 4])][t,ρ]-180 (t^4) ρ Derivative[1,1][(Subscript[f, 4])][t,ρ]-210 t ρ Derivative[2,0][ξ][t,ρ]+315 (t^4) ρ Derivative[2,0][(Subscript[b, 4])][t,ρ])))+(1/(3780 (r^7) (t^7) (ρ^2)))(17016 (ρ^2)+105840 (t^5) (ρ^2) ((ξ[t,ρ])^5)+47880 (t^6) (ρ^2) ((ξ[t,ρ])^6)+1152 (t^4) (ρ^2) (Subscript[a, 4])[t,ρ]+84 (t^4) (ρ^2) ((ξ[t,ρ])^4) (1931-1350 (t^2) Derivative[1,0][ξ][t,ρ])+6 (t^3) ρ ((ξ[t,ρ])^3) (30529 ρ+210 (t^2) (-60 (t^2) ρ (Subscript[b, 4])[t,ρ]+t Derivative[0,1][ξ][t,ρ]+ρ (t Derivative[0,2][ξ][t,ρ]-129 Derivative[1,0][ξ][t,ρ]+21 t Derivative[2,0][ξ][t,ρ])))-6 (t^2) ρ ((ξ[t,ρ])^2) (-24274 ρ-1260 (t^4) ρ (Subscript[a, 4])[t,ρ]-1260 (t^5) (Subscript[f, 4])[t,ρ]+(t^2) (103 t Derivative[0,1][ξ][t,ρ]-420 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+2520 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]-277 t ρ Derivative[0,2][ξ][t,ρ]+24861 ρ Derivative[1,0][ξ][t,ρ]-8820 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)-9450 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]-165 (t^2) Derivative[1,1][ξ][t,ρ]+15 (t^2) ρ Derivative[1,2][ξ][t,ρ]-5250 t ρ Derivative[2,0][ξ][t,ρ]+315 (t^2) ρ Derivative[3,0][ξ][t,ρ]))+3 ξ[t,ρ] (24206 t (ρ^2)+2352 (t^5) (ρ^2) (Subscript[a, 4])[t,ρ]+(t^3) (84 (t^2) (-15 (t^2)+16 (ρ^2)) (Subscript[b, 4])[t,ρ]+84 (15 (t^4)+16 (t^2) (ρ^2)) (Subscript[c, 4])[t,ρ]+ρ (1764 (t^3) (Subscript[f, 4])[t,ρ]+688 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+147 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-1050 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]+840 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]+1764 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+304 t ρ Derivative[0,2][ξ][t,ρ]+336 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-210 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-27462 ρ Derivative[1,0][ξ][t,ρ]-60 (t^3) ρ Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+12180 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)-630 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]+330 (t^2) Derivative[1,1][ξ][t,ρ]+12 t Derivative[0,1][ξ][t,ρ] (-68+5 (t^2) (11 Derivative[1,0][ξ][t,ρ]-2 ρ Derivative[1,1][ξ][t,ρ]))+2520 (t^4) ρ Derivative[1,1][(Subscript[f, 4])][t,ρ]-30 (t^2) ρ Derivative[1,2][ξ][t,ρ]+6510 t ρ Derivative[2,0][ξ][t,ρ]-3780 (t^3) ρ Derivative[1,0][ξ][t,ρ] Derivative[2,0][ξ][t,ρ]-4410 (t^4) ρ Derivative[2,0][(Subscript[b, 4])][t,ρ]-630 (t^2) ρ Derivative[3,0][ξ][t,ρ])))-(t^2) (-12 (35 (t^4)+54 (t^2) (ρ^2)) (Subscript[b, 4])[t,ρ]+12 (35 (t^4)-54 (t^2) (ρ^2)) (Subscript[c, 4])[t,ρ]+120 (t^3) ρ (Subscript[f, 4])[t,ρ]+1286 t ρ Derivative[0,1][ξ][t,ρ]-492 (t^2) (ρ^2) ((Derivative[0,1][ξ][t,ρ])^2)+93 (t^4) ρ Derivative[0,1][(Subscript[a, 4])][t,ρ]-378 (t^4) ρ Derivative[0,1][(Subscript[b, 4])][t,ρ]-66 (t^4) ρ Derivative[0,1][(Subscript[c, 4])][t,ρ]-930 (t^3) (ρ^2) Derivative[0,1][(Subscript[f, 4])][t,ρ]-194 t (ρ^2) Derivative[0,2][ξ][t,ρ]+78 (t^4) (ρ^2) Derivative[0,2][(Subscript[a, 4])][t,ρ]+96 (t^4) (ρ^2) Derivative[0,2][(Subscript[b, 4])][t,ρ]+60 (t^4) (ρ^2) Derivative[0,2][(Subscript[c, 4])][t,ρ]+21162 (ρ^2) Derivative[1,0][ξ][t,ρ]-990 (t^3) ρ Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+180 (t^4) (ρ^2) ((Derivative[0,1][ξ][t,ρ])^2) Derivative[1,0][ξ][t,ρ]+90 (t^3) (ρ^2) Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]-10710 (t^2) (ρ^2) ((Derivative[1,0][ξ][t,ρ])^2)+3780 (t^4) (ρ^2) ((Derivative[1,0][ξ][t,ρ])^3)-90 (t^5) Derivative[1,0][(Subscript[b, 4])][t,ρ]+351 (t^3) (ρ^2) Derivative[1,0][(Subscript[b, 4])][t,ρ]+90 (t^5) Derivative[1,0][(Subscript[c, 4])][t,ρ]+36 (t^3) (ρ^2) Derivative[1,0][(Subscript[c, 4])][t,ρ]-330 (t^2) ρ Derivative[1,1][ξ][t,ρ]+180 (t^3) (ρ^2) Derivative[0,1][ξ][t,ρ] Derivative[1,1][ξ][t,ρ]-90 (t^5) ρ Derivative[1,1][(Subscript[a, 4])][t,ρ]+135 (t^5) ρ Derivative[1,1][(Subscript[b, 4])][t,ρ]+360 (t^5) ρ Derivative[1,1][(Subscript[c, 4])][t,ρ]+210 (t^4) (ρ^2) Derivative[1,1][(Subscript[f, 4])][t,ρ]+30 (t^2) (ρ^2) Derivative[1,2][ξ][t,ρ]+180 (t^5) (ρ^2) Derivative[1,2][(Subscript[a, 4])][t,ρ]-45 (t^5) (ρ^2) Derivative[1,2][(Subscript[b, 4])][t,ρ]-5460 t (ρ^2) Derivative[2,0][ξ][t,ρ]+5670 (t^3) (ρ^2) Derivative[1,0][ξ][t,ρ] Derivative[2,0][ξ][t,ρ]-288 (t^4) (ρ^2) Derivative[2,0][(Subscript[a, 4])][t,ρ]-315 (t^4) (ρ^2) Derivative[2,0][(Subscript[b, 4])][t,ρ]-390 (t^5) ρ Derivative[2,0][(Subscript[f, 4])][t,ρ]+780 (t^5) (ρ^2) Derivative[2,1][(Subscript[f, 4])][t,ρ]+630 (t^2) (ρ^2) Derivative[3,0][ξ][t,ρ]-945 (t^5) (ρ^2) Derivative[3,0][(Subscript[b, 4])][t,ρ]))/.r->1/z;
Cseries[z_,t_,ρ_]=-(2/(3 r t))+Log[(1/(t^(2/3)))]+((1+2 t ξ[t,ρ])/(3 (r^2) (t^2)))-((2 (1+3 t ξ[t,ρ]+3 (t^2) ((ξ[t,ρ])^2)))/(9 (r^3) (t^3)))+(((Subscript[c, 4])[t,ρ])/(r^4))+(((8/(15 (t^5)))+((2 ξ[t,ρ])/(t^4))+((10 ((ξ[t,ρ])^2))/(3 (t^3)))+((10 ((ξ[t,ρ])^3))/(3 (t^2)))+((2 ((ξ[t,ρ])^4))/t)+((2 (Subscript[a, 4])[t,ρ])/(15 t))-4 ξ[t,ρ] (Subscript[c, 4])[t,ρ]-((4 (Subscript[f, 4])[t,ρ])/(15 ρ))+(2/15) Derivative[0,1][(Subscript[f, 4])][t,ρ]-((2 Derivative[1,0][ξ][t,ρ])/(3 (t^3)))-((2 ξ[t,ρ] Derivative[1,0][ξ][t,ρ])/(t^2))-((2 ((ξ[t,ρ])^2) Derivative[1,0][ξ][t,ρ])/t)+Derivative[1,0][(Subscript[c, 4])][t,ρ])/(r^5))+(1/(540 (r^6) (t^6) ρ))(-960 ρ+t (-6300 (t^3) ρ ((ξ[t,ρ])^4)-3240 (t^4) ρ ((ξ[t,ρ])^5)+90 (t^2) ρ ((ξ[t,ρ])^3) (-91+60 (t^2) Derivative[1,0][ξ][t,ρ])+30 ((ξ[t,ρ])^2) (180 (t^5) ρ (Subscript[c, 4])[t,ρ]-(t^4) Derivative[0,1][ξ][t,ρ]+t ρ (-240+(t^2) (5 t Derivative[0,2][ξ][t,ρ]+219 Derivative[1,0][ξ][t,ρ]-21 t Derivative[2,0][ξ][t,ρ])))-30 ξ[t,ρ] (129 ρ+12 (t^4) ρ (Subscript[a, 4])[t,ρ]-24 (t^5) (Subscript[f, 4])[t,ρ]+(t^2) (t Derivative[0,1][ξ][t,ρ]-10 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+ρ (-141 Derivative[1,0][ξ][t,ρ]+t (-5 Derivative[0,2][ξ][t,ρ]+6 t (2 t Derivative[0,1][(Subscript[f, 4])][t,ρ]+7 ((Derivative[1,0][ξ][t,ρ])^2)+15 t Derivative[1,0][(Subscript[c, 4])][t,ρ])+21 Derivative[2,0][ξ][t,ρ]))))+t (-72 (t^2) ρ (Subscript[a, 4])[t,ρ]-30 (t^3) (Subscript[f, 4])[t,ρ]-10 t Derivative[0,1][ξ][t,ρ]+150 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)-24 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]+30 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]-15 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]+60 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+50 t ρ Derivative[0,2][ξ][t,ρ]+12 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-30 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-45 (t^4) ρ Derivative[0,2][(Subscript[c, 4])][t,ρ]+1230 ρ Derivative[1,0][ξ][t,ρ]-630 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+72 (t^3) ρ Derivative[1,0][(Subscript[a, 4])][t,ρ]+45 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]-180 (t^4) Derivative[1,0][(Subscript[f, 4])][t,ρ]+90 (t^4) ρ Derivative[1,1][(Subscript[f, 4])][t,ρ]-210 t ρ Derivative[2,0][ξ][t,ρ]+315 (t^4) ρ Derivative[2,0][(Subscript[c, 4])][t,ρ])))+(1/(7560 (r^7) (t^7) (ρ^2)))(34032 (ρ^2)+211680 (t^5) (ρ^2) ((ξ[t,ρ])^5)+95760 (t^6) (ρ^2) ((ξ[t,ρ])^6)+2304 (t^4) (ρ^2) (Subscript[a, 4])[t,ρ]+168 (t^4) (ρ^2) ((ξ[t,ρ])^4) (1931-1350 (t^2) Derivative[1,0][ξ][t,ρ])-12 (t^3) ρ ((ξ[t,ρ])^3) (-30529 ρ+210 (t^2) (60 (t^2) ρ (Subscript[c, 4])[t,ρ]+11 t Derivative[0,1][ξ][t,ρ]+ρ (5 t Derivative[0,2][ξ][t,ρ]+129 Derivative[1,0][ξ][t,ρ]-21 t Derivative[2,0][ξ][t,ρ])))+12 (t^2) ρ ((ξ[t,ρ])^2) (24274 ρ+1260 (t^4) ρ (Subscript[a, 4])[t,ρ]-2520 (t^5) (Subscript[f, 4])[t,ρ]+(t^2) (-2173 t Derivative[0,1][ξ][t,ρ]-2100 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+1260 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]-1223 t ρ Derivative[0,2][ξ][t,ρ]-24861 ρ Derivative[1,0][ξ][t,ρ]+8820 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+9450 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]-15 (t^2) Derivative[1,1][ξ][t,ρ]+165 (t^2) ρ Derivative[1,2][ξ][t,ρ]+5250 t ρ Derivative[2,0][ξ][t,ρ]-315 (t^2) ρ Derivative[3,0][ξ][t,ρ]))-12 ξ[t,ρ] (-12103 t (ρ^2)-1176 (t^5) (ρ^2) (Subscript[a, 4])[t,ρ]+(t^3) (-84 (15 (t^4)+8 (t^2) (ρ^2)) (Subscript[b, 4])[t,ρ]+84 (15 (t^4)-8 (t^2) (ρ^2)) (Subscript[c, 4])[t,ρ]+ρ (378 (t^3) (Subscript[f, 4])[t,ρ]+1396 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+147 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-1050 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]-105 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]-252 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+688 t ρ Derivative[0,2][ξ][t,ρ]+84 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-210 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-315 (t^4) ρ Derivative[0,2][(Subscript[c, 4])][t,ρ]+13731 ρ Derivative[1,0][ξ][t,ρ]-330 (t^3) ρ Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]-6090 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+315 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]+15 (t^2) Derivative[1,1][ξ][t,ρ]+6 t Derivative[0,1][ξ][t,ρ] (103+5 (t^2) (Derivative[1,0][ξ][t,ρ]-22 ρ Derivative[1,1][ξ][t,ρ]))+630 (t^4) ρ Derivative[1,1][(Subscript[f, 4])][t,ρ]-165 (t^2) ρ Derivative[1,2][ξ][t,ρ]-3255 t ρ Derivative[2,0][ξ][t,ρ]+1890 (t^3) ρ Derivative[1,0][ξ][t,ρ] Derivative[2,0][ξ][t,ρ]+2205 (t^4) ρ Derivative[2,0][(Subscript[c, 4])][t,ρ]+315 (t^2) ρ Derivative[3,0][ξ][t,ρ])))+(t^2) ((-420 (t^4)+1296 (t^2) (ρ^2)) (Subscript[b, 4])[t,ρ]+12 (35 (t^4)+108 (t^2) (ρ^2)) (Subscript[c, 4])[t,ρ]+2280 (t^3) ρ (Subscript[f, 4])[t,ρ]+668 t ρ Derivative[0,1][ξ][t,ρ]-4056 (t^2) (ρ^2) ((Derivative[0,1][ξ][t,ρ])^2)-51 (t^4) ρ Derivative[0,1][(Subscript[a, 4])][t,ρ]-684 (t^4) ρ Derivative[0,1][(Subscript[b, 4])][t,ρ]-48 (t^4) ρ Derivative[0,1][(Subscript[c, 4])][t,ρ]+600 (t^3) (ρ^2) Derivative[0,1][(Subscript[f, 4])][t,ρ]-2012 t (ρ^2) Derivative[0,2][ξ][t,ρ]+24 (t^4) (ρ^2) Derivative[0,2][(Subscript[a, 4])][t,ρ]-12 (t^4) (ρ^2) Derivative[0,2][(Subscript[b, 4])][t,ρ]+60 (t^4) (ρ^2) Derivative[0,2][(Subscript[c, 4])][t,ρ]-42324 (ρ^2) Derivative[1,0][ξ][t,ρ]-180 (t^3) ρ Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+3960 (t^4) (ρ^2) ((Derivative[0,1][ξ][t,ρ])^2) Derivative[1,0][ξ][t,ρ]+1980 (t^3) (ρ^2) Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+21420 (t^2) (ρ^2) ((Derivative[1,0][ξ][t,ρ])^2)-7560 (t^4) (ρ^2) ((Derivative[1,0][ξ][t,ρ])^3)-180 (t^5) Derivative[1,0][(Subscript[b, 4])][t,ρ]-72 (t^3) (ρ^2) Derivative[1,0][(Subscript[b, 4])][t,ρ]+180 (t^5) Derivative[1,0][(Subscript[c, 4])][t,ρ]-702 (t^3) (ρ^2) Derivative[1,0][(Subscript[c, 4])][t,ρ]-60 (t^2) ρ Derivative[1,1][ξ][t,ρ]+3960 (t^3) (ρ^2) Derivative[0,1][ξ][t,ρ] Derivative[1,1][ξ][t,ρ]-360 (t^5) ρ Derivative[1,1][(Subscript[a, 4])][t,ρ]+360 (t^5) ρ Derivative[1,1][(Subscript[b, 4])][t,ρ]-270 (t^5) ρ Derivative[1,1][(Subscript[c, 4])][t,ρ]+840 (t^4) (ρ^2) Derivative[1,1][(Subscript[f, 4])][t,ρ]+660 (t^2) (ρ^2) Derivative[1,2][ξ][t,ρ]+180 (t^5) (ρ^2) Derivative[1,2][(Subscript[a, 4])][t,ρ]-360 (t^5) (ρ^2) Derivative[1,2][(Subscript[b, 4])][t,ρ]-630 (t^5) (ρ^2) Derivative[1,2][(Subscript[c, 4])][t,ρ]+10920 t (ρ^2) Derivative[2,0][ξ][t,ρ]-11340 (t^3) (ρ^2) Derivative[1,0][ξ][t,ρ] Derivative[2,0][ξ][t,ρ]+576 (t^4) (ρ^2) Derivative[2,0][(Subscript[a, 4])][t,ρ]+630 (t^4) (ρ^2) Derivative[2,0][(Subscript[c, 4])][t,ρ]-1560 (t^5) ρ Derivative[2,0][(Subscript[f, 4])][t,ρ]+780 (t^5) (ρ^2) Derivative[2,1][(Subscript[f, 4])][t,ρ]-1260 (t^2) (ρ^2) Derivative[3,0][ξ][t,ρ]+1890 (t^5) (ρ^2) Derivative[3,0][(Subscript[c, 4])][t,ρ]))/.r->1/z;
Sseries[z_,t_,ρ_]=-(1/(9 r (t^(5/3))))+r (t^(1/3))+((1+3 t ξ[t,ρ])/(3 (t^(2/3))))+((5+9 t ξ[t,ρ])/(81 (r^2) (t^(8/3))))+((-10-3 t ξ[t,ρ] (10+9 t ξ[t,ρ]))/(243 (r^3) (t^(11/3))))+((440 ρ+9 t (200 ρ ξ[t,ρ]+300 t ρ ((ξ[t,ρ])^2)+180 (t^2) ρ ((ξ[t,ρ])^3)-27 (t^3) (4 ρ (Subscript[a, 4])[t,ρ]+t (4 (Subscript[f, 4])[t,ρ]+4 ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+3 ρ Derivative[1,0][(Subscript[a, 4])][t,ρ]))))/(14580 (r^4) (t^(14/3)) ρ))+(1/(131220 (r^5) (t^(17/3)) ρ))(2752 ρ+9 t (-2952 (t^2) ρ ((ξ[t,ρ])^3)-1620 (t^3) ρ ((ξ[t,ρ])^4)+72 ((ξ[t,ρ])^2) (-23 t ρ+9 (t^3) (t Derivative[0,1][ξ][t,ρ]+t ρ Derivative[0,2][ξ][t,ρ]-3 ρ Derivative[1,0][ξ][t,ρ]))+4 ξ[t,ρ] (46 ρ+81 (t^2) (2 t (Derivative[0,1][ξ][t,ρ]+2 t (3 ρ (Subscript[a, 4])[t,ρ]+ρ ((Derivative[0,1][ξ][t,ρ])^2)+3 t ((Subscript[f, 4])[t,ρ]+ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]))+ρ Derivative[0,2][ξ][t,ρ])-6 ρ Derivative[1,0][ξ][t,ρ]+9 (t^3) ρ Derivative[1,0][(Subscript[a, 4])][t,ρ]))+27 t (24 (t^2) ρ (Subscript[a, 4])[t,ρ]+12 (t^3) (Subscript[f, 4])[t,ρ]+8 t Derivative[0,1][ξ][t,ρ]+24 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+3 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-24 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]+12 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]+12 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+8 t ρ Derivative[0,2][ξ][t,ρ]+3 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-12 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-24 ρ Derivative[1,0][ξ][t,ρ]-9 (t^3) ρ Derivative[1,0][(Subscript[a, 4])][t,ρ]+18 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]+18 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]-9 (t^4) ρ Derivative[2,0][(Subscript[a, 4])][t,ρ])))+(1/(5511240 (r^6) (t^(20/3)) (ρ^2)))(-1342588 (ρ^2)-6055560 (t^4) (ρ^2) ((ξ[t,ρ])^4)-4286520 (t^5) (ρ^2) ((ξ[t,ρ])^5)-236196 (t^4) (ρ^2) (Subscript[a, 4])[t,ρ]-648 (t^3) ρ ((ξ[t,ρ])^3) (8809 ρ+1890 (t^2) (t Derivative[0,1][ξ][t,ρ]+t ρ Derivative[0,2][ξ][t,ρ]-3 ρ Derivative[1,0][ξ][t,ρ]))+216 (t^2) ρ ((ξ[t,ρ])^2) (-23863 ρ+135 (t^2) (126 (t^2) ρ (Subscript[b, 4])[t,ρ]+126 (t^2) ρ (Subscript[c, 4])[t,ρ]-62 t Derivative[0,1][ξ][t,ρ]-84 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)-47 t ρ Derivative[0,2][ξ][t,ρ]+195 ρ Derivative[1,0][ξ][t,ρ]+9 (t^2) Derivative[1,1][ξ][t,ρ]+9 (t^2) ρ Derivative[1,2][ξ][t,ρ]-27 t ρ Derivative[2,0][ξ][t,ρ]))-6 t ρ ξ[t,ρ] (652952 ρ+1215 (t^2) (252 (t^2) ρ (Subscript[a, 4])[t,ρ]+84 (t^2) ρ (Subscript[b, 4])[t,ρ]+84 (t^2) ρ (Subscript[c, 4])[t,ρ]+168 (t^3) (Subscript[f, 4])[t,ρ]+172 t Derivative[0,1][ξ][t,ρ]+208 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2)+21 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]-168 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]+84 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]+168 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+112 t ρ Derivative[0,2][ξ][t,ρ]+21 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-84 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-660 ρ Derivative[1,0][ξ][t,ρ]-72 (t^3) Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]-72 (t^3) ρ Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+216 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+126 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]+126 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]-36 (t^2) Derivative[1,1][ξ][t,ρ]-144 (t^3) ρ Derivative[0,1][ξ][t,ρ] Derivative[1,1][ξ][t,ρ]-36 (t^2) ρ Derivative[1,2][ξ][t,ρ]+108 t ρ Derivative[2,0][ξ][t,ρ]-63 (t^4) ρ Derivative[2,0][(Subscript[a, 4])][t,ρ]))+243 (t^2) (12 (75 (t^4)-11 (t^2) (ρ^2)) (Subscript[b, 4])[t,ρ]-12 (75 (t^4)+11 (t^2) (ρ^2)) (Subscript[c, 4])[t,ρ]+ρ (-1644 (t^3) (Subscript[f, 4])[t,ρ]-219 (t^4) Derivative[0,1][(Subscript[a, 4])][t,ρ]+690 (t^4) Derivative[0,1][(Subscript[b, 4])][t,ρ]-30 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]-744 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]-920 t ρ Derivative[0,2][ξ][t,ρ]+6 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-150 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]-90 (t^4) ρ Derivative[0,2][(Subscript[c, 4])][t,ρ]+7080 ρ Derivative[1,0][ξ][t,ρ]+1080 (t^3) ρ Derivative[0,2][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]-3240 (t^2) ρ ((Derivative[1,0][ξ][t,ρ])^2)+240 (t^2) ρ ((Derivative[0,1][ξ][t,ρ])^2) (-7+9 (t^2) Derivative[1,0][ξ][t,ρ])-450 (t^3) ρ Derivative[1,0][(Subscript[b, 4])][t,ρ]-450 (t^3) ρ Derivative[1,0][(Subscript[c, 4])][t,ρ]+360 (t^2) Derivative[1,1][ξ][t,ρ]+40 t Derivative[0,1][ξ][t,ρ] (-38+27 (t^2) (Derivative[1,0][ξ][t,ρ]+2 ρ Derivative[1,1][ξ][t,ρ]))+135 (t^5) Derivative[1,1][(Subscript[a, 4])][t,ρ]-1080 (t^5) Derivative[1,1][(Subscript[b, 4])][t,ρ]+540 (t^5) Derivative[1,1][(Subscript[c, 4])][t,ρ]+900 (t^4) ρ Derivative[1,1][(Subscript[f, 4])][t,ρ]+360 (t^2) ρ Derivative[1,2][ξ][t,ρ]+135 (t^5) ρ Derivative[1,2][(Subscript[a, 4])][t,ρ]-540 (t^5) ρ Derivative[1,2][(Subscript[b, 4])][t,ρ]-1080 t ρ Derivative[2,0][ξ][t,ρ]+495 (t^4) ρ Derivative[2,0][(Subscript[a, 4])][t,ρ]+810 (t^4) ρ Derivative[2,0][(Subscript[b, 4])][t,ρ]+810 (t^4) ρ Derivative[2,0][(Subscript[c, 4])][t,ρ]+720 (t^5) Derivative[2,0][(Subscript[f, 4])][t,ρ]+720 (t^5) ρ Derivative[2,1][(Subscript[f, 4])][t,ρ]+135 (t^5) ρ Derivative[3,0][(Subscript[a, 4])][t,ρ])))/.r->1/z;
Fseries[z_,t_,ρ_]=(((Subscript[f, 4])[t,ρ])/(r^2))+Derivative[0,1][ξ][t,ρ]+((-10 ξ[t,ρ] (Subscript[f, 4])[t,ρ]+Derivative[0,1][(Subscript[a, 4])][t,ρ]+4 Derivative[1,0][(Subscript[f, 4])][t,ρ])/(5 (r^3)))+(1/(r^4))(3 ((ξ[t,ρ])^2) (Subscript[f, 4])[t,ρ]+((Derivative[0,1][ξ][t,ρ])/(18 (t^4)))-((((ξ[t,ρ])^2) Derivative[0,1][ξ][t,ρ])/(6 (t^2)))+((Derivative[0,1][(Subscript[a, 4])][t,ρ])/(12 t))-(3/5) ξ[t,ρ] Derivative[0,1][(Subscript[a, 4])][t,ρ]+((Derivative[0,1][(Subscript[b, 4])][t,ρ])/(12 t))+((Derivative[0,1][(Subscript[c, 4])][t,ρ])/(12 t))-((Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ])/(6 (t^2)))-((ξ[t,ρ] Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ])/(3 t))+((Derivative[1,0][(Subscript[b, 4])][t,ρ])/(12 ρ))-((Derivative[1,0][(Subscript[c, 4])][t,ρ])/(12 ρ))-(12/5) ξ[t,ρ] Derivative[1,0][(Subscript[f, 4])][t,ρ]-((Derivative[1,1][ξ][t,ρ])/(18 (t^3)))-((ξ[t,ρ] Derivative[1,1][ξ][t,ρ])/(6 (t^2)))-((((ξ[t,ρ])^2) Derivative[1,1][ξ][t,ρ])/(6 t))+(1/6) Derivative[1,1][(Subscript[a, 4])][t,ρ]+(1/12) Derivative[1,1][(Subscript[b, 4])][t,ρ]+(1/3) Derivative[2,0][(Subscript[f, 4])][t,ρ])-(1/(1260 (r^5) (t^5) (ρ^2)))(12096 (t^4) (ρ^2) ((ξ[t,ρ])^4) Derivative[0,1][ξ][t,ρ]+168 (t^3) (ρ^2) ((ξ[t,ρ])^3) (30 (t^2) (Subscript[f, 4])[t,ρ]+67 Derivative[0,1][ξ][t,ρ]-5 t Derivative[1,1][ξ][t,ρ])-12 (t^2) ((ξ[t,ρ])^2) (504 (t^3) ρ (Subscript[b, 4])[t,ρ]-504 (t^3) ρ (Subscript[c, 4])[t,ρ]-504 (t^2) (ρ^2) (Subscript[f, 4])[t,ρ]-3 (t^2) Derivative[0,1][ξ][t,ρ]-312 (ρ^2) Derivative[0,1][ξ][t,ρ]+504 (t^3) (ρ^2) Derivative[0,1][(Subscript[b, 4])][t,ρ]+3 (t^2) ρ Derivative[0,2][ξ][t,ρ]+3 (t^2) (ρ^2) Derivative[0,3][ξ][t,ρ]+140 (t^2) (ρ^2) Derivative[0,1][ξ][t,ρ] Derivative[1,0][ξ][t,ρ]+76 t (ρ^2) Derivative[1,1][ξ][t,ρ]-15 (t^2) (ρ^2) Derivative[2,1][ξ][t,ρ])+4 t ξ[t,ρ] (-18 (t^3) ρ ((Derivative[0,1][ξ][t,ρ])^2)+Derivative[0,1][ξ][t,ρ] (9 (t^2)+106 (ρ^2)-6 (t^2) (ρ^2) (9 t Derivative[0,2][ξ][t,ρ]+41 Derivative[1,0][ξ][t,ρ]-15 t Derivative[2,0][ξ][t,ρ]))+t ρ (105 (t^2) ρ Derivative[0,1][(Subscript[a, 4])][t,ρ]+105 (t^2) ρ Derivative[0,1][(Subscript[b, 4])][t,ρ]+105 (t^2) ρ Derivative[0,1][(Subscript[c, 4])][t,ρ]-9 t Derivative[0,2][ξ][t,ρ]-9 t ρ Derivative[0,3][ξ][t,ρ]+105 (t^3) Derivative[1,0][(Subscript[b, 4])][t,ρ]-105 (t^3) Derivative[1,0][(Subscript[c, 4])][t,ρ]-178 ρ Derivative[1,1][ξ][t,ρ]+180 (t^2) ρ Derivative[1,0][ξ][t,ρ] Derivative[1,1][ξ][t,ρ]+210 (t^3) ρ Derivative[1,1][(Subscript[a, 4])][t,ρ]+105 (t^3) ρ Derivative[1,1][(Subscript[b, 4])][t,ρ]+420 (t^3) ρ Derivative[2,0][(Subscript[f, 4])][t,ρ]+45 t ρ Derivative[2,1][ξ][t,ρ]))+3 (12 (t^4) (Subscript[f, 4])[t,ρ]-12 (t^3) ρ ((Derivative[0,1][ξ][t,ρ])^2)-24 (t^4) (ρ^2) ((Derivative[0,1][ξ][t,ρ])^3)+4 Derivative[0,1][ξ][t,ρ] ((t^2)+26 (ρ^2)+3 (t^2) (ρ^2) (-3 t Derivative[0,2][ξ][t,ρ]+2 Derivative[1,0][ξ][t,ρ] (-6+5 (t^2) Derivative[1,0][ξ][t,ρ])+5 t Derivative[2,0][ξ][t,ρ]))-t (3 ((t^4)-10 (t^2) (ρ^2)) Derivative[0,1][(Subscript[a, 4])][t,ρ]+6 (2 (t^4)-5 (t^2) (ρ^2)) Derivative[0,1][(Subscript[b, 4])][t,ρ]-6 (t^4) Derivative[0,1][(Subscript[c, 4])][t,ρ]-30 (t^2) (ρ^2) Derivative[0,1][(Subscript[c, 4])][t,ρ]+12 (t^3) ρ Derivative[0,1][(Subscript[f, 4])][t,ρ]+4 t ρ Derivative[0,2][ξ][t,ρ]-3 (t^4) ρ Derivative[0,2][(Subscript[a, 4])][t,ρ]-12 (t^4) ρ Derivative[0,2][(Subscript[b, 4])][t,ρ]+6 (t^4) ρ Derivative[0,2][(Subscript[c, 4])][t,ρ]+12 (t^3) (ρ^2) Derivative[0,2][(Subscript[f, 4])][t,ρ]+4 t (ρ^2) Derivative[0,3][ξ][t,ρ]-3 (t^4) (ρ^2) Derivative[0,3][(Subscript[a, 4])][t,ρ]-6 (t^4) (ρ^2) Derivative[0,3][(Subscript[b, 4])][t,ρ]+88 (ρ^2) Derivative[1,1][ξ][t,ρ]-120 (t^2) (ρ^2) Derivative[1,0][ξ][t,ρ] Derivative[1,1][ξ][t,ρ]+33 (t^3) (ρ^2) Derivative[1,1][(Subscript[a, 4])][t,ρ]+24 (t^3) (ρ^2) Derivative[1,1][(Subscript[b, 4])][t,ρ]+24 (t^3) (ρ^2) Derivative[1,1][(Subscript[c, 4])][t,ρ]+30 (t^4) ρ Derivative[2,0][(Subscript[b, 4])][t,ρ]-30 (t^4) ρ Derivative[2,0][(Subscript[c, 4])][t,ρ]-20 t (ρ^2) Derivative[2,1][ξ][t,ρ]+33 (t^4) (ρ^2) Derivative[2,1][(Subscript[a, 4])][t,ρ]+30 (t^4) (ρ^2) Derivative[2,1][(Subscript[b, 4])][t,ρ]+40 (t^4) (ρ^2) Derivative[3,0][(Subscript[f, 4])][t,ρ])))/.r->1/z;

Sdotseries[z_,t_,ρ_]=Collect[Normal[Series[D[Sseries[z,t,ρ],t]+(1/2)Aseries[z,t,ρ]D[Sseries[1/r,t,ρ],r]/.r->1/z,{z,0,4}]],z,Simplify];
Bdotseries[z_,t_,ρ_]=Collect[Normal[Series[D[Bseries[z,t,ρ],t]+(1/2)Aseries[z,t,ρ]D[Bseries[1/r,t,ρ],r]/.r->1/z,{z,0,4}]],z,Simplify];
Cdotseries[z_,t_,ρ_]=Collect[Normal[Series[D[Cseries[z,t,ρ],t]+(1/2)Aseries[z,t,ρ]D[Cseries[1/r,t,ρ],r]/.r->1/z,{z,0,4}]],z,Simplify];
Fdotseries[z_,t_,ρ_]=Collect[Normal[Series[D[Fseries[z,t,ρ],t]+(1/2)Aseries[z,t,ρ]D[Fseries[1/r,t,ρ],r]/.r->1/z,{z,0,4}]],z,Simplify];

replf4=(Subscript[f, 4])->((Subscript[fs, 4])[#1,#2]#2&);

{a4p[t_,ρ_],f4sp[t_,ρ_]}={-((4 (-ρ-4 t ρ ξ[t,ρ]-6 (t^2) ρ ((ξ[t,ρ])^2)-4 (t^3) ρ ((ξ[t,ρ])^3)+3 (t^4) (t (Subscript[f, 4])[t,ρ]+ρ ((Subscript[a, 4])[t,ρ]+(Subscript[b, 4])[t,ρ]+(Subscript[c, 4])[t,ρ]+t Derivative[0,1][(Subscript[f, 4])][t,ρ]))))/(9 (t^5) ρ)),-((-12 (t^3) (Subscript[b, 4])[t,ρ]+12 (t^3) (Subscript[c, 4])[t,ρ]+ρ (12 (t^2) (Subscript[f, 4])[t,ρ]+8 (1+3 t ξ[t,ρ] (1+t ξ[t,ρ])) Derivative[0,1][ξ][t,ρ]+3 (t^3) (Derivative[0,1][(Subscript[a, 4])][t,ρ]-4 Derivative[0,1][(Subscript[b, 4])][t,ρ])))/(12 (t^3) ρ))(1/ρ)}/.replf4//FullSimplify

Series[f4sp[t,ρ],{ρ,0,1}];

bi[z_,ρ_]:=(E^(-(ρ^2))) (z^4);
ci[z_,ρ_]:=(E^(-(ρ^2))) (z^4);
a4i[ρ_]=-10 (E^(-(ρ^2)));
f4i[ρ_]=0;
bi2[z_,ρ_]:=-(1/4)(E^(-(ρ^2))) (z^4)+(11/30)(E^(-(ρ^2))) (z^5)(E^(-(z^2)))/ti;
ci2[z_,ρ_]:=-(1/4)(E^(-(ρ^2))) (z^4)+(11/30)(E^(-(ρ^2))) (z^5)(E^(-(z^2)))/ti;

bi3[z_,ρ_]:=-(1/4) (z^4)+(11/30) (z^5)(E^(-(z^2)))/ti;
ci3[z_,ρ_]:=-(1/4) (z^4)+(11/30) (z^5)(E^(-(z^2)))/ti;
a4i3[ρ_]=-10 ;

Initialise numerics

Initialises various replacements rules to construct formulas, grid and numerical functions etc. Not everything is relevant. Particularly important: 
changeN(n): changes number of gridpoints
InitNum(t0, profile, Binitial, a4initial), initialises numerics (also called by run)

It may be useful to realise that the code internally uses truncated series, i.e. it works with Bfin(z,t) instead of B(z,t), where B(z,t)=((∑_(i=0))^(powerB-1))(Subscript[b, i])(z^i)+(z^powerB)Bfin(z,t), and analogously for all the other functions. This is the relevance of setting powerB, powerS etc.

Also it is good to realise that to use analytic functions of z and t, one should first call ClearRepl, since in the numerics z and t get numeric values.

Useful shortcuts

machinep=True;
If[machinep==False,
$MinPrecision=50;,
$MinPrecision=0;]

N2[exp_]:=Block[{$MinPrecision=5},N[exp 1.,5]];
N3[exp_]:=Block[{$MinPrecision=5},N[Max[Abs[exp]] 1.,5]];
M2[exp_]:=Block[{$MinPrecision=5},MatrixForm[N[exp 1.,5]]];
list3D[funct_]:=Table[{z[[i]],ρ[[k]],funct[[k,i]]},{i,nz},{k,2,nρ}]
(*pl[exp_]:=ListPlot3D[Flatten[list3D[exp],1],PlotRange->All,AxesLabel->{"z","ρ"}]*)
pl[exp_]:=(fint=Interpolation[Flatten[list3D[exp],1]];
Plot3D[fint[zp,ρp],{zp,0,z[[1]]},{ρp,0,ρ[[2]]},PlotPoints->30,PlotRange->All,AxesLabel->{"z","ρ"}])
pllist[list_,func_]:=Show[pl2[list,{func[0],func[40]}],Plot[func[u],{u,0,4L},PlotRange->All],ImageSize->350];
plnums[deriv_:1,n_:10]:=(nums[[#,deriv]]//pl)&/@Range[n];
findRange[list_,minrange_:0]:=(
max=Max[list];
If[minrange=!=0,max=Max[{max,minrange[[1]],minrange[[2]]}]];
min=Min[list];
If[minrange=!=0,min=Min[{min,minrange[[1]],minrange[[2]]}]];
{min-(max-min)/10,max+(max-min)/10}
);
pl2[exp_,minrange_:0]:=ListPlot[Thread[{{1}~Join~((ρ[[2;;-1]]/L)/(1+ρ[[2;;-1]]/L)),exp}],ImageSize->300,PlotRange->{{0,1.01},Evaluate[findRange[exp[[Length[Select[ρ,#>L 10&]];;-1]],minrange]]}]
pl2[exp_,minrange_:0]:=ListPlot[Thread[{{ρ[[2]] 1.05}~Join~ρ[[2;;-1]],exp}],ImageSize->300,PlotRange->{{0,ρ[[2]]1.1},Evaluate[findRange[exp[[1;;-1]],minrange]]}]
pllist[list_,func_]:=Show[pl2[list],Plot[func[u],{u,0,ρ[[2]]1.1},PlotRange->All],ImageSize->350]
plz[exp_]:=ListPlot[Thread[{z,exp}],PlotRange->All]
plρ[exp_]:=pl2[exp];
tρ[exp_]:=Thread[{ρ,exp}];
Dr[func_,times_]:=D[func[1/r,t,ρ],{r,times}]/.r->1/z
(*replori={(Subscript[f, 4])->(0&),(Subscript[b, 4])->(Subscript[c, 4]),Derivative[a_,b_][f_][t,0]:> If[OddQ[b],0,Derivative[a,b][f][t,0]]};*)
replori[exp_]:=exp/.{Derivative[a_,b_][f_][t,0]:>If[OddQ[b],0,Derivative[a,b][f][t,0]]}/.{(Subscript[b, 4])->(Subscript[c, 4]),(Subscript[f, 4])->(0&)}
Simplify2[exp_]:=Simplify[exp/.ρ->(ρ^3),Assumptions->{ρ>0,t>0}]/.ρ->(ρ^(1/3))
FullSimplify2[exp_]:=FullSimplify[exp/.ρ->(ρ^3),Assumptions->{ρ>0,t>0}]/.ρ->(ρ^(1/3))

ClearRepl:=Module[{},
Unset[{zp,tp,x,g,z,ρ,t,a4,b4,c4,i,fs4,fs4dρ,b4dρ,a4dρ,c4dρ,ksi,ksidρ,ksidρ2,ksidt,ksidρdt,ksidtdρ,a4dt,fs4dt,ksidρ2dt}];
Unset[{Bf,Bfd,Ff,Ffd,Sf,Sfd,Sdotf,Bdotf,Fnum,Bnumdρ,Bnumdz,Bnumdρdz,Af,nums,numsb,addksiT}];
Unset[{ρpos,zpos,rp,rpp,trial,trialdρ,trialdρ2,tabCheb,D2Specρ,DSpecρ,DSpec}];
];
cl:=ClearRepl;
cln:=(Clear[nz];Clear[nρ];ClearRepl;)
cln

New grid, init and update

points=7;
order=4;
Λ=0;
changeN[newnz_,newnρ_]:=Module[{},
If[newnz>27||newnρ>27||True,
machinep=False;
$MinPrecision=60;
$MaxPrecision=2 $MinPrecision;,
machinep=True;
$MinPrecision=0];
nz=newnz;
nρ=2 newnρ-1;
If[machinep,
zCheb=N[Table[Cos[((r π)/(nz-1))],{r,0,nz-1}],30];
ρCheb=N[Table[Cos[((r π)/(nρ-1))],{r,0,nρ-1}],30];
ρFour=N[Table[((π ii)/(nρ-1)),{ii,0,nρ-1}]];,
zCheb=N[Table[Cos[((r π)/(nz-1))],{r,0,nz-1}],$MinPrecision];
ρCheb=N[Table[Cos[((r π)/(nρ-1))],{r,0,nρ-1}],$MinPrecision];
ρFour=N[Table[((π ii)/(nρ-1)),{ii,0,nρ-1}],$MinPrecision]];
Cheb=Table[ChebyshevT[i-1,zCheb[[k]]],{i,nz},{k,nz}];
Chebρ=Table[ChebyshevT[i-1,ρCheb[[k]]],{i,nρ},{k,nρ}];
dSpec=Table[If[i==nz&&j==nz,-((2(((nz-1))^2)+1)/6),If[i==1&&j==1,((2(((nz-1))^2)+1)/6),If[i==j,((-zCheb[[i]])/(2(1-((zCheb[[i]])^2)))),((((-1))^(i+j))/(zCheb[[i]]-zCheb[[j]]))((If[i==1||i==nz,2,1])/(If[j==1||j==nz,2,1]))]]],{i,nz},{j,nz}];
Chebρ=Table[ChebyshevT[i-1,ρCheb[[k]]],{i,nρ},{k,nρ}];
dSpecρ=Table[If[i==nρ&&j==nρ,-((2(((nρ-1))^2)+1)/6),If[i==1&&j==1,((2(((nρ-1))^2)+1)/6),If[i==j,((-ρCheb[[i]])/(2(1-((ρCheb[[i]])^2)))),((((-1))^(i+j))/(ρCheb[[i]]-ρCheb[[j]]))((If[i==1||i==nρ,2,1])/(If[j==1||j==nρ,2,1]))]]],{i,nρ},{j,nρ}];
d2Specρ=dSpecρ.dSpecρ;
dFour=Table[If[i==1||i==nρ,0,If[i==j,Cot[ρFour[[i]]]/2,(((-1))^j)((Sin[ρFour[[i]]]Cos[(nρ-1)ρFour[[i]]])/((Subscript[pcardρ, j]) (Cos[ρFour[[i]]]-Cos[ρFour[[j]]])))]],{i, nρ},{j, nρ}];
dFourSin=Table[If[j==1||j==nρ,0,If[i==j,-Cot[ρFour[[i]]]/2,(((-1))^(j+i+1))((Sin[ρFour[[j]]])/((Cos[ρFour[[i]]]-Cos[ρFour[[j]]])))]],{i, nρ},{j, nρ}];
dFour2=Table[If[i==j,If[machinep,-(1+2 (((nρ/2))^2))/6.,N[-(1+2 (((nρ/2))^2))/6,$MinPrecision]],(((-1))^(i-j+1))(1/(2((Sin[(ρFour[[i]]-ρFour[[j]])/2])^2)))],{i,nρ},{j,nρ}];
dFour2=dFourSin.dFour;
nρ=newnρ;
];

ClearRepl;
add0[m_]:=Table[If[k≤Length[m[[1]]],m[[i,k]],0],{i,Length[m]},{k,Length[m[[1]]]+1}];
newGrid[lowerz_,upperρ_]:=Module[{},
ϵz=lowerz;
ϵρ=upperρ;
ax=(ϵz-Λ) /2;
bx=(ϵz+Λ)/2;
axρ=(ϵρ-Λ) /2;
bxρ=(ϵρ+Λ)/2;
z=ax zCheb+bx;
zs=z[[1;;-2]];
ρ=axρ ρCheb+bxρ;
ρ=(E^((((axρ ρCheb+bxρ))^3)+(axρ ρCheb+bxρ)))-1;
ρ=((ϵρ)/π)ρFour//Reverse;
L=ϵρ;
ρ=L((ρCheb)/(Sqrt[1-((ρCheb)^2)]))[[1;;nρ]];
ρ=L(((((ρCheb)/(((1-((ρCheb)^2)))^(1/powercoord)))))^1)[[1;;nρ]];
ρ[[1]]=∞;
ρs=ρ[[1;;-2]];
(*zplist=Transpose[Table[Table[If[j==n,0,replz[[i,1]]/.z->zx[[j]]],{j,n}],{i,Length[replz]}]];*)
DSpec=dSpec/ax;
D2Spec = DSpec.DSpec;
(*DSpecρ=(dSpecρ/axρ);
DSpecρ=((E^(-(((axρ ρCheb+bxρ))^2)-(axρ ρCheb+bxρ)))/((3(((axρ ρCheb+bxρ))^2)+1)axρ)) dSpecρ;
D2Specρ = DSpecρ.DSpecρ;
DSpecρ=(π/(ϵρ))Table[dFour[[nρ-i+1,nρ-k+1]],{i,nρ},{k,nρ}];
D2Specρ=((π^2)/((ϵρ)^2))Table[dFour2[[nρ-i+1,nρ-k+1]],{i,nρ},{k,nρ}];
DSpecρ=DiagonalMatrix[((((1-((ρCheb)^2)))^(3/2))/L)[[1;;nρ]]].dSpecρ[[1;;nρ,1;;nρ]]+DiagonalMatrix[((((1-((ρCheb)^2)))^(3/2))/L)[[1;;nρ]]].add0[dSpecρ[[1;;nρ,-1;;nρ+1;;-1]]];
D2Specρ=DiagonalMatrix[((((1-((ρCheb)^2)))^3)/(L^2))[[1;;nρ]]].d2Specρ[[1;;nρ,1;;nρ]]+DiagonalMatrix[((((1-((ρCheb)^2)))^3)/(L^2))[[1;;nρ]]].add0[d2Specρ[[1;;nρ,-1;;nρ+1;;-1]]]-3( DiagonalMatrix[(((((1-((ρCheb)^2)))^2)ρCheb)/(L^2))[[1;;nρ]]].dSpecρ[[1;;nρ,1;;nρ]]+DiagonalMatrix[(((((1-((ρCheb)^2)))^2)ρCheb)/(L^2))[[1;;nρ]]].add0[dSpecρ[[1;;nρ,-1;;nρ+1;;-1]]]);*)
DSpecρ=DiagonalMatrix[(( powercoord (((1-((ρCheb)^2)))^(1+(1/powercoord))))/(L powercoord+2 L ((ρCheb)^2)-L powercoord ((ρCheb)^2)))[[1;;nρ]]].dSpecρ[[1;;nρ,1;;nρ]]+DiagonalMatrix[(( powercoord (((1-((ρCheb)^2)))^(1+(1/powercoord))))/(L powercoord+2 L ((ρCheb)^2)-L powercoord ((ρCheb)^2)))[[1;;nρ]]].add0[dSpecρ[[1;;nρ,-1;;nρ+1;;-1]]];
D2Specρ=DiagonalMatrix[-(( (powercoord^2) (((1-((ρCheb)^2)))^(1+(2/powercoord))) (-1+((ρCheb)^2)))/((L^2) (((-powercoord+(-2+powercoord) ((ρCheb)^2)))^2)))[[1;;nρ]]].d2Specρ[[1;;nρ,1;;nρ]]+DiagonalMatrix[-(( (powercoord^2) (((1-((ρCheb)^2)))^(1+(2/powercoord))) (-1+((ρCheb)^2)))/((L^2) (((-powercoord+(-2+powercoord) ((ρCheb)^2)))^2)))[[1;;nρ]]].add0[d2Specρ[[1;;nρ,-1;;nρ+1;;-1]]]+( DiagonalMatrix[-(( powercoord (((1-((ρCheb)^2)))^(1+(2/powercoord))) (-6 powercoord ρCheb+2 (-2+powercoord) ((ρCheb)^3)))/((L^2) (((-powercoord+(-2+powercoord) ((ρCheb)^2)))^3)))[[1;;nρ]]].dSpecρ[[1;;nρ,1;;nρ]]+DiagonalMatrix[-(( powercoord (((1-((ρCheb)^2)))^(1+(2/powercoord))) (-6 powercoord ρCheb+2 (-2+powercoord) ((ρCheb)^3)))/((L^2) (((-powercoord+(-2+powercoord) ((ρCheb)^2)))^3)))[[1;;nρ]]].add0[dSpecρ[[1;;nρ,-1;;nρ+1;;-1]]]);
(*constant=(a/der2)/.Solve[(DSpecρ.{4a,2a,a}~Join~Table[0,{nρ-3}])[[1]]==der2,a][[1]];*)
{N2[z],N2[ρ]}
];

{dρ,dρ2,dz,dz2,dzdρ,dρdz}={2,3,4,5,6,6};
repl[funct_,n_]:={funct[z,t,ρ]->nums[[n,1,i]],Derivative[0,0,1][funct][z,t,ρ]->nums[[n,2,i]],Derivative[0,0,2][funct][z,t,ρ]->nums[[n,3,i]],Derivative[1,0,0][funct][z,t,ρ]->nums[[n,4,i]],Derivative[2,0,0][funct][z,t,ρ]->nums[[n,5,i]],Derivative[1,0,1][funct][z,t,ρ]->nums[[n,6,i]],Derivative[0,0,4][funct][z,t,0]->numsb[[n,1]],Derivative[1,0,2][funct][z,t,0]->numsb[[n,2]]}

powers={2,2,1,1,0,1,1,1};

funcs={B,C,S,F,Sdot,Cdot,Bdot,A,Sdotdt,Fdt};
{nB,nC,nS,nF,nSdot,nCdot,nBdot,nA,nSdotdt,nFdt}=Range[10];
powers={2,2,1,0,0,1,1,0};
powers={3,3,3,1,1,2,2,1};

updateCoef:=(
t=time;
If[ksidtfunc=!= (0&),
ksi=ksifunc[t,#]&/@ρ;
ksidt=ksidtfunc[t,#]&/@ρ;];
If[powers[[nB]]==3,
b4=Transpose[nums[[nB,1]]][[nz]],
b4=Transpose[nums[[nB,dz]]][[nz]]];
If[powers[[nC]]==3,
c4=Transpose[nums[[nC,1]]][[nz]],
c4=Transpose[nums[[nC,dz]]][[nz]]];
b4dρ=DSpecρ.b4;
c4dρ=DSpecρ.c4;
a4dρ=DSpecρ.a4;
fs4dρ=DSpecρ.fs4;
ksidρ=DSpecρ.ksi;
ksidρ2=D2Specρ.ksi;
ksidρdt=DSpecρ.ksidt;
ksidρ2dt=D2Specρ.ksidt;)

update[n_]:=(
nums[[n,dz]]=(DSpec.#1&)/@nums[[n,1]];
nums[[n,dz2]]=(D2Spec.#1&)/@nums[[n,1]];
nums[[n,dρ]]=DSpecρ.nums[[n,1]];
nums[[n,dρ2]]=D2Specρ.nums[[n,1]];
nums[[n,dzdρ]]=DSpecρ.nums[[n,dz]];
numsb[[n,1]]=(D2Specρ.nums[[n,dρ2]])[[nρ]];
numsb[[n,2]]=DSpec.nums[[n,dρ2,nρ]];)
SetAttributes[update,Listable]

getnum[func_,ksiini_]:=(
boundz[ρp_,tp_]=SeriesCoefficient[Series[func[ztmp,tp,ρp],{ztmp,0,0}],0];
boundρ[zp_,tp_]=SeriesCoefficient[Series[func[zp,tp,ρtmp],{ρtmp,0,0}],0]/.Log[ρtmp]->0;
boundzρ=SeriesCoefficient[Series[func[ztmp,t,ρtmp],{ρtmp,0,0},{ztmp,0,0}],{0,0}];
(*((func[zs,t,xtmp]/.xtmp->#)~Join~{boundz[#,t]}&/@ρs)~Join~
{boundρ[zs,t]~Join~{boundzρ}}/.ξ->ksiini*)
((func[zs,t,xtmp]/.xtmp->#)~Join~{boundz[#,t]}&/@ρ)/.ξ->ksiini
)

BAdS[z_,t_,ρ_]:=-(2/3)(Log[(t+z)](*+Log[ρ]*));
CAdS[z_,t_,ρ_]:=-(2/3)Log[(t+z)(*/(ρ^2)*)];
init[Tini_,profB_,profC_,f4ini_,a4ini_,ksiini_,zsize_:1,ρsize_:1]:=(
ti=Tini;
time=ti;
t=time;
newGrid[zsize,ρsize];
nums=Table[,{nFdt},{dρdz}];
numsb=Table[,{nFdt},{2}];
Binitot[z_,t_,ρ_]=(BAdS[(z/(1+z ξ[t,ρ])),t,ρ]+profB[(z/(1+z ξ[t,ρ])),ρ]-((Subscript[subtract, B])[z,t,ρ]/.Log[t ρ]->Log[t]))/(z^(powers[[nB]]+1));
Cinitot[z_,t_,ρ_]=(CAdS[(z/(1+z ξ[t,ρ])),t,ρ]+profC[(z/(1+z ξ[t,ρ])),ρ]-((Subscript[subtract, C])[z,t,ρ]/.Log[((ρ^(4/3))/(t^(2/3)))]->-(2/3)Log[t]))/(z^(powers[[nC]]+1));
nums[[nB,1]]=getnum[Binitot,ksiini];
nums[[nC,1]]=getnum[Cinitot,ksiini];
fs4=(f4ini[#]/#&/@ρs)~Join~{If[Series[f4ini[ρp],{ρp,0,1}]==0,0,SeriesCoefficient[Series[f4ini[ρp],{ρp,0,1}],1]]};
a4=a4ini/@ρ;
ksifunc=ksiini;
ksidtfunc=Derivative[1,0][ksiini];
ksi=ksifunc[t,#]&/@ρ;
ksidt=ksidtfunc[t,#]&/@ρ;
update[{nB,nC}];
updateCoef;
b4ini=Interpolation[Thread[{ρ,b4}][[2;;-1]]];
c4ini=Interpolation[Thread[{ρ,c4}][[2;;-1]]];)

Do[(Subscript[subtract, funcs[[i]]])[z_,t_,ρ_]=Normal[Series[ToExpression[ToString[funcs[[i]]]<>"series[z,t,ρ]"],{z,0,powers[[i]]}]]//FullSimplify,{i,Length[powers]}]

replfunc=Table[funcs[[i]]->(Evaluate[(Subscript[subtract, Evaluate[funcs[[i]]]])[#1,#2,#3]+(If[i==nBdot,If[i==nBdot,(((#3)^2)/(1+((#3)^2))),1] ToExpression[ToString[funcs[[i]]]<>"f[#1,#2,#3]"]+ToExpression[ToString[funcs[[nCdot]]]<>"f[#1,#2,#3]"],ToExpression[ToString[funcs[[i]]]<>"f[#1,#2,#3]"]])((#1)^(powers[[i]]+1))If[i==nS||i==nSdot,1,If[i==nF,((#3)/(1+((#3)^2))),1]]]&),{i,Length[powers]}];
replfunc=Table[replfunc[[i]],{i,Length[replfunc]}]

Off[Part::pspec]
Off[Part::partd]

replcoef={ρ->ρ[[i]],(Subscript[a, 4])[t,ρ]->a4[[i]],(Subscript[b, 4])[t,ρ]->b4[[i]],(Subscript[c, 4])[t,ρ]->c4[[i]],(Subscript[fs, 4])[t,ρ]->fs4[[i]],Derivative[0,1][(Subscript[fs, 4])][t,ρ]->fs4dρ[[i]],Derivative[0,1][(Subscript[b, 4])][t,ρ]->b4dρ[[i]],Derivative[0,1][(Subscript[c, 4])][t,ρ]->c4dρ[[i]],Derivative[0,1][(Subscript[a, 4])][t,ρ]->a4dρ[[i]],Derivative[1,0][(Subscript[fs, 4])][t,ρ]->fs4dt[[i]],Derivative[1,0][(Subscript[b, 4])][t,ρ]->b4dt[[i]],Derivative[1,0][(Subscript[c, 4])][t,ρ]->c4dt[[i]],Derivative[1,0][(Subscript[a, 4])][t,ρ]->a4dt[[i]],ξ[t,ρ]->ksi[[i]],Derivative[1,0][ξ][t,ρ]->ksidt[[i]],Derivative[0,1][ξ][t,ρ]->ksidρ[[i]],Derivative[0,2][ξ][t,ρ]->ksidρ2[[i]],Derivative[1,0][ξ][t,ρ]->ksidt[[i]],Derivative[1,1][ξ][t,ρ]->ksidρdt[[i]],Derivative[1,2][ξ][t,ρ]->ksidρ2dt[[i]],Derivative[0,4][ξ][t,0]->(D2Specρ.ksidρ2)[[nρ]], Derivative[2,0,1][f_][z,t,0]->0,Derivative[1,0,2][Bf][z,t,0]->DSpec.nums[[nB,dρ2,-1]],Derivative[1,0,2][Sf][z,t,0]->DSpec.nums[[nS,dρ2,-1]],Derivative[1,0,2][Cf][z,t,0]->DSpec.nums[[nC,dρ2,-1]]};

replnum=Flatten[Table[repl[ToExpression[ToString[funcs[[ii]]]<>"f"],ii],{ii,Length[funcs]}]]
replnumtoint={nums[[nn_,der_,ρi_]]->  interpolateCheb[nn,der,ρi,zpos]}
replnumtointρz=nums[[nC_,der_,ρi_]]->interpolateChebFour[nC,der,ρpos,zpos];
replintksi={ksi[[i]]->intksi[ρpos],ksidρ[[i]]->intksidρ[ρpos],ksidρ2[[i]]->intksidρ2[ρpos]};
replnumb=replnum/.z->0/.nums[[a_,b_,i]]-> nums[[a,b,i,nz]];
replnumbρ=replnum~Join~replcoef[[2;;-1]]/.ρ->0/.nums[[a_,b_,i]]-> nums[[a,b,nρ]]/.f_[[i]]->f[[nρ]];

Interpolation, filtering of high modes and changing the grid

(Subscript[pcard, i_]):=If[i==nz||i==1,2,1];
(Subscript[pcardρ, i_]):=If[i==nρ||i==1,2,1];
chebnz=-1;
(Subscript[CardCheb, i_])[xlow_]:=(x=SetPrecision[xlow,Max[$MinPrecision,100]];If[x==zCheb[[i]],1,(((-1))^i)(((1-(x^2)))/((Subscript[pcard, i])(((nz-1))^2)(x-zCheb[[i]])))(nz-1)If[chebnz==nz&&chebx==x,cheb,chebnz=nz;chebx=x;cheb=ChebyshevU[nz-2,x]]]);
(Subscript[Card, i_])[x_,ax_,bx_]:=(Subscript[CardCheb, i])[(x-bx)/ax];

interpolateCheb[nC_,der_,ypos_Integer,zpos_]:=Sum[(Subscript[Card, j])[zpos,ax,bx]nums[[nC,der,ypos,j]],{j,nz}]
interpolateFour[nC_,der_,ypos_,zpos_Integer]:=(fourtmp=Re[Fourier[nums[[nC,der,1;;2 ny,zpos]]]];(((∑_i^ny)fourtmp[[i]] Cos[((2 π  (((2 ny)/ymax)ypos) (i-1))/(2 ny))]+fourtmp[[ny+1]] Cos[(((2 ny)/ymax)ypos) π]+(∑_(i=ny+2)^(2 ny))fourtmp[[i]] Cos[-((2 π  (((2 ny)/ymax)ypos) (2 ny-i+1))/(2 ny))])/(Sqrt[2 ny])))
interpolateFourOdd[nC_,der_,ypos_,zpos_Integer]:=(fourtmp=Im[Fourier[nums[[nC,der,1;;2 ny,zpos]]]];(((∑_i^ny)fourtmp[[i]] Sin[((2 π  (((2 ny)/ymax)ypos) (i-1))/(2 ny))]+fourtmp[[ny+1]] Sin[(((2 ny)/ymax)ypos) π]+(∑_(i=ny+2)^(2 ny))fourtmp[[i]] Sin[-((2 π  (((2 ny)/ymax)ypos) (2 ny-i+1))/(2 ny))])/(Sqrt[2 ny])))
orderint=9;
interpolateFourFast[nC_,der_,ypos_,zpos_Integer]:=Interpolation[Thread[{y,nums[[nC,der,1;;2 ny,zpos]]}],InterpolationOrder->9][ypos]
intksi[ypos_Real]:=Interpolation[Thread[{y,ksi}],InterpolationOrder->9][ypos]
intksidy[ypos_Real]:=Interpolation[Thread[{y,ksidy}],InterpolationOrder->9][ypos]
intksidy2[ypos_Real]:=Interpolation[Thread[{y,ksidy2}],InterpolationOrder->9][ypos]
interpolateChebFour[nC_,der_,ypos_,zpos_Real]:=(
indy=Round[1+2ny(ypos/ymax)];
yindmin=Max[1,Min[indy-Round[orderint/2],2 ny-orderint-1]];
yindmax=Max[orderint+1,Min[indy-Round[orderint/2]+orderint,2 ny]];
Interpolation[Thread[{y[[yindmin;;yindmax]],(interpolateCheb[nC,der,#,zpos]&/@Range[yindmin,yindmax])}],InterpolationOrder->orderint][ypos])

Nmodes=2;
filterFour[list_]:=(
fourtmp=Re[Fourier[list]];
list-Table[(((∑_(i=ny-Nmodes+1)^ny)fourtmp[[i]] Cos[((2 π  (((2 ny)/ymax)y[[k]]) (i-1))/(2 ny))]+fourtmp[[ny+1]] Cos[((((2 ny) y[[k]]) π)/ymax)]+(∑_(i=ny+2)^(ny+Nmodes+1))fourtmp[[i]] Cos[-((2 π  (((2 ny)/ymax)y[[k]]) (2 ny-i+1))/(2 ny))])/(Sqrt[2 ny])),{k,2 ny}]);
filterFourEven[list_]:=filterFour[list~Join~list[[-2;;-ny;;-1]]][[1;;ny+1]]
filterFourOdd[list_]:=filterFour[list~Join~(-list[[-2;;-ny;;-1]])][[1;;ny+1]]

nfilterz=0;
nfilterρ=0;
nfilterρ2=0;

filterCheb[list_]:=(
coef=Table[(2/(nz-1))Sum[(1/(Subscript[pcard, k]))(list[[k]])Cheb[[i,k]],{k,nz}],{i,nz-nfilterz+1,nz}];
list-Table[Sum[Cheb[[nz-i+1,i]] coef[[-i2]],{i2,nfilterz}],{i, nz}]

);
filterChebρ[list_]:=(
(*coef=Table[(2/(nρ-1))Sum[(1/(Subscript[pcardρ, k]))(list[[k]])Chebρ[[i,k]],{k,nρ}],{i,nρ-5,nρ}];
list-Table[2 Chebρ[[nρ-2,i]] coef[[1]]/3+3Chebρ[[nρ-1,i]] coef[[2]]/4+5 Chebρ[[nρ,i]] coef[[3]]/6,{i,nρ}]*)
coefeven=Table[(2/((2 nρ-1)-1))(1/(Subscript[pcardρ, i]))(Sum[(2/(Subscript[pcardρ, k]))list[[k]]Chebρ[[2 i-1,k]],{k,nρ}]),{i,nρ-nfilterρ+1,nρ}];
list-Table[Sum[Chebρ[[2 nρ-2 i2+1,i]] coefeven[[-i2]],{i2,nfilterρ}],{i, nρ}]
);
nfilterρ2:=0;
filterChebρ2[list_]:=(
(*coef=Table[(2/(nρ-1))Sum[(1/(Subscript[pcardρ, k]))(list[[k]])Chebρ[[i,k]],{k,nρ}],{i,nρ-5,nρ}];
list-Table[2 Chebρ[[nρ-2,i]] coef[[1]]/3+3Chebρ[[nρ-1,i]] coef[[2]]/4+5 Chebρ[[nρ,i]] coef[[3]]/6,{i,nρ}]*)
coefeven=Table[(2/((2 nρ-1)-1))(1/(Subscript[pcardρ, i]))(Sum[(2/(Subscript[pcardρ, k]))list[[k]]Chebρ[[2 i-1,k]],{k,nρ}]),{i,nρ-nfilterρ2+1,nρ}];
list-Table[Sum[Chebρ[[2 nρ-2 i2+1,i]] coefeven[[-i2]],{i2,nfilterρ2}],{i, nρ}]9/10
);
filterBoundDer[list_]:=(
(*der=(DSpecρ.list)[[1]];
list-({4 constant der,2constant der,constant der}~Join~Table[0,{nρ-3}])*)
list)

(*Table[(Subscript[Card, i])[z[[1]]1.0005,ax,bx],{i,nz}];
Sum[%[[i]],{i,nz}];*)

trans[list_,addksi_,power_,sub_]:=Table[(CTmatrix[[i]].list[[i]])(((1+addksiT[[i]] z))^(-power))+sub,{i,nρ}];
changeGrid[addksi_]:=(
addksiT=addksi;
CTmatrix=Table[(Subscript[Card, j])[((z[[i]])/(1+z[[i]] addksi[[k]])),ax,bx],{k,nρ},{i,nz},{j,nz}];
nums[[nB,1]]=trans[nums[[nB,1]],addksi,powers[[nB]]+1,subB];
nums[[nC,1]]=trans[nums[[nC,1]],addksi,powers[[nC]]+1,subC];
update[{nB,nC}];
Bdtlist=trans[#,addksi,powers[[nB]]+1,subBdt]&/@Bdtlist;
Cdtlist=trans[#,addksi,powers[[nC]]+1,subBdt]&/@Cdtlist;
Sdotlist=trans[#,addksi,powers[[nSdot]]+1,subSdot]&/@Sdotlist;
Flist=trans[#,addksi,powers[[nF]]+1,subF]&/@Flist;
ksi+=addksi;
updateCoef;
);

subB=(1/(z^(powers[[nB]]+1)))((Subscript[subtract, B])[(z/(1+z addksiT[[i]])),t,ρ]-((Subscript[subtract, B])[z,t,ρ]/.ξ[t,ρ]->ξ[t,ρ]+addksiT[[i]]))/.replcoef//Simplify
subBdt=D[subB,t]//Simplify
subC=(1/(z^(powers[[nC]]+1)))((Subscript[subtract, C])[(z/(1+z addksiT[[i]])),t,ρ]-((Subscript[subtract, C])[z,t,ρ]/.ξ[t,ρ]->ξ[t,ρ]+addksiT[[i]]))/.replcoef//Simplify
subSdot=(1/(z^(powers[[nSdot]]+1)))((Subscript[subtract, Sdot])[(z/(1+z addksiT[[i]])),t,ρ]-((Subscript[subtract, Sdot])[z,t,ρ]/.ξ[t,ρ]->ξ[t,ρ]+addksiT[[i]]))/.replcoef//Simplify
subF=(1/(z^(powers[[nF]]+1)))((Subscript[subtract, F])[(z/(1+z addksiT[[i]])),t,ρ]-((Subscript[subtract, F])[z,t,ρ]/.ξ[t,ρ]->ξ[t,ρ]+addksiT[[i]]))/.replcoef//Simplify

Full equations

Compute equations

ClearRepl
compute=False;
session="fullrhoandFb";

replrtoz={Derivative[ddr_,ddt_,ddρ_][f_][r,t,ρ]->Hold[Dr[Derivative[0,ddt,ddρ][f],ddr]],r->z}

D[Derivative[0,0,1][f][r,t,y]-F[r,t,y] Derivative[1,0,0][f][r,t,y],y]-F[r,t,y]D[(Derivative[0,0,1][f][r,t,y]-F[r,t,y] Derivative[1,0,0][f][r,t,y]),r]//Simplify

D[Derivative[0,0,1][f][r,t,y]-F[r,t,y] Derivative[1,0,0][f][r,t,y],r]//FullSimplify

D[Derivative[0,1,0][f][r,t,y]+(1/2) A[r,t,y] Derivative[1,0,0][f][r,t,y],y]-F[r,t,y]D[Derivative[0,1,0][f][r,t,y]+(1/2) A[r,t,y] Derivative[1,0,0][f][r,t,y],r]//FullSimplify
D[Derivative[0,1,0][f][r,t,y]+(1/2) A[r,t,y] Derivative[1,0,0][f][r,t,y],r]//FullSimplify

repldρ={Derivative[0,0,2][f_][r,t,ρ]->Derivative[0,0,2][f][r,t,ρ]-Derivative[0,0,1][F][r,t,ρ] Derivative[1,0,0][f][r,t,ρ]+F[r,t,ρ] (Derivative[1,0,0][f][r,t,ρ] Derivative[1,0,0][F][r,t,ρ]-2 Derivative[1,0,1][f][r,t,ρ]+F[r,t,ρ] Derivative[2,0,0][f][r,t,ρ]),
Derivative[0,1,1][f_][r,t,ρ]->(1/2) (2 Derivative[0,1,1][f][r,t,ρ]+Derivative[0,0,1][A][r,t,ρ] Derivative[1,0,0][f][r,t,ρ]+A[r,t,ρ] Derivative[1,0,1][f][r,t,ρ]-F[r,t,ρ] (Derivative[1,0,0][A][r,t,ρ] Derivative[1,0,0][f][r,t,ρ]+2 Derivative[1,1,0][f][r,t,ρ]+A[r,t,ρ] Derivative[2,0,0][f][r,t,ρ])) ,
Derivative[1,1,0][f_][r,t,ρ]->(1/2) (Derivative[1,0,0][A][r,t,ρ] Derivative[1,0,0][f][r,t,ρ]+2 Derivative[1,1,0][f][r,t,ρ]+A[r,t,ρ] Derivative[2,0,0][f][r,t,ρ]),
Derivative[1,0,1][f_][r,t,ρ]->-Derivative[1,0,0][f][r,t,ρ] Derivative[1,0,0][F][r,t,ρ]+Derivative[1,0,1][f][r,t,ρ]-F[r,t,ρ] Derivative[2,0,0][f][r,t,ρ],
Derivative[0,0,1][f_][r,t,ρ]->Derivative[0,0,1][f][r,t,ρ]-F[r,t,ρ] Derivative[1,0,0][f][r,t,ρ],
Derivative[0,1,0][F][r,t,ρ]->Derivative[0,1,0][F][r,t,ρ]+(1/2) A[r,t,ρ] Derivative[1,0,0][F][r,t,ρ],
Derivative[0,1,0][A][r,t,ρ]->Derivative[0,1,0][A][r,t,ρ]+(1/2) A[r,t,ρ] Derivative[1,0,0][A][r,t,ρ]};
repldot={Derivative[ddr_,1,ddρ_][B][r_,t,ρ]->Derivative[ddr,0,ddρ][Bdot][r,t,ρ],Derivative[ddr_,1,ddρ_][C][r_,t,ρ]->Derivative[ddr,0,ddρ][Cdot][r,t,ρ],Derivative[ddr_,1,ddρ_][S][r_,t,ρ]->Derivative[ddr,0,ddρ][Sdot][r,t,ρ](*,Derivative[ddr_,1,ddρ_][F][r_,t,ρ]->Derivative[ddr,0,ddρ][Fdot][r,t,ρ]*),
Derivative[0,2,0][S][r,t,ρ]->Derivative[0,1,0][Sdot][r,t,ρ]+(1/2) A[r,t,ρ] Derivative[1,0,0][Sdot][r,t,ρ]};

EE={Derivative[2,0,0][S][r,t,ρ]->-(1/6) S[r,t,ρ] (((Derivative[1,0,0][B][r,t,ρ])^2)+Derivative[1,0,0][B][r,t,ρ] Derivative[1,0,0][C][r,t,ρ]+((Derivative[1,0,0][C][r,t,ρ])^2)),Derivative[2,0,0][F][r,t,ρ]->(1/(2 ρ ((S[r,t,ρ])^2)))(-8 ρ Derivative[0,0,1][S][r,t,ρ] Derivative[1,0,0][S][r,t,ρ]+((S[r,t,ρ])^2) ((-2+2 ρ Derivative[0,0,1][B][r,t,ρ]+ρ Derivative[0,0,1][C][r,t,ρ]) Derivative[1,0,0][B][r,t,ρ]+(2+ρ Derivative[0,0,1][B][r,t,ρ]+2 ρ Derivative[0,0,1][C][r,t,ρ]) Derivative[1,0,0][C][r,t,ρ]-2 ρ Derivative[1,0,1][B][r,t,ρ])+2 ρ S[r,t,ρ] (-3 Derivative[0,0,1][S][r,t,ρ] Derivative[1,0,0][B][r,t,ρ]+3 Derivative[1,0,0][F][r,t,ρ] Derivative[1,0,0][S][r,t,ρ]+4 Derivative[1,0,1][S][r,t,ρ])),Derivative[1,1,0][S][r,t,ρ]->(1/(12 ρ ((S[r,t,ρ])^3)))(E^(-B[r,t,ρ])) (24 (E^(B[r,t,ρ])) ρ ((S[r,t,ρ])^4)+4 ρ ((Derivative[0,0,1][S][r,t,ρ])^2)+2 S[r,t,ρ] (-4 ρ Derivative[0,0,2][S][r,t,ρ]+Derivative[0,0,1][S][r,t,ρ] (-4+4 ρ Derivative[0,0,1][B][r,t,ρ]+ρ Derivative[1,0,0][F][r,t,ρ]))-((S[r,t,ρ])^2) (2 ρ ((Derivative[0,0,1][B][r,t,ρ])^2)+2 Derivative[0,0,1][C][r,t,ρ]+ρ ((Derivative[0,0,1][C][r,t,ρ])^2)-2 Derivative[1,0,0][F][r,t,ρ]+Derivative[0,0,1][B][r,t,ρ] (-4+ρ Derivative[0,0,1][C][r,t,ρ]+2 ρ Derivative[1,0,0][F][r,t,ρ])-ρ (2 Derivative[0,0,2][B][r,t,ρ]+((Derivative[1,0,0][F][r,t,ρ])^2)-24 (E^(B[r,t,ρ])) Derivative[0,1,0][S][r,t,ρ] Derivative[1,0,0][S][r,t,ρ]+2 Derivative[1,0,1][F][r,t,ρ]))),Derivative[1,1,0][C][r,t,ρ]->(1/(6 ρ ((S[r,t,ρ])^4)))(E^(-B[r,t,ρ])) (-4 ρ ((Derivative[0,0,1][S][r,t,ρ])^2)+S[r,t,ρ] (2 ρ Derivative[0,0,2][S][r,t,ρ]-Derivative[0,0,1][S][r,t,ρ] (4+2 ρ Derivative[0,0,1][B][r,t,ρ]+3 ρ Derivative[0,0,1][C][r,t,ρ]-4 ρ Derivative[1,0,0][F][r,t,ρ]))-9 (E^(B[r,t,ρ])) ρ ((S[r,t,ρ])^3) (Derivative[0,1,0][S][r,t,ρ] Derivative[1,0,0][C][r,t,ρ]+Derivative[0,1,0][C][r,t,ρ] Derivative[1,0,0][S][r,t,ρ])+((S[r,t,ρ])^2) (2 ρ ((Derivative[0,0,1][B][r,t,ρ])^2)+ρ ((Derivative[0,0,1][C][r,t,ρ])^2)+4 Derivative[1,0,0][F][r,t,ρ]+2 Derivative[0,0,1][B][r,t,ρ] (1+2 ρ Derivative[0,0,1][C][r,t,ρ]+ρ Derivative[1,0,0][F][r,t,ρ])+Derivative[0,0,1][C][r,t,ρ] (-1+3 ρ Derivative[1,0,0][F][r,t,ρ])-ρ (2 Derivative[0,0,2][B][r,t,ρ]+3 Derivative[0,0,2][C][r,t,ρ]+((Derivative[1,0,0][F][r,t,ρ])^2)+2 Derivative[1,0,1][F][r,t,ρ]))),Derivative[1,1,0][B][r,t,ρ]->(1/(6 ρ ((S[r,t,ρ])^4)))(E^(-B[r,t,ρ])) (8 ρ ((Derivative[0,0,1][S][r,t,ρ])^2)+S[r,t,ρ] (-4 ρ Derivative[0,0,2][S][r,t,ρ]+Derivative[0,0,1][S][r,t,ρ] (2+ρ Derivative[0,0,1][B][r,t,ρ]-8 ρ Derivative[1,0,0][F][r,t,ρ]))-9 (E^(B[r,t,ρ])) ρ ((S[r,t,ρ])^3) (Derivative[0,1,0][S][r,t,ρ] Derivative[1,0,0][B][r,t,ρ]+Derivative[0,1,0][B][r,t,ρ] Derivative[1,0,0][S][r,t,ρ])-((S[r,t,ρ])^2) (ρ ((Derivative[0,0,1][B][r,t,ρ])^2)+4 Derivative[0,0,1][C][r,t,ρ]+2 ρ ((Derivative[0,0,1][C][r,t,ρ])^2)+2 Derivative[1,0,0][F][r,t,ρ]+Derivative[0,0,1][B][r,t,ρ] (1+2 ρ Derivative[0,0,1][C][r,t,ρ]+ρ Derivative[1,0,0][F][r,t,ρ])-ρ (Derivative[0,0,2][B][r,t,ρ]+2 ((Derivative[1,0,0][F][r,t,ρ])^2)+4 Derivative[1,0,1][F][r,t,ρ]))),Derivative[2,0,0][A][r,t,ρ]->(1/(2 ρ ((S[r,t,ρ])^4)))(E^(-B[r,t,ρ])) (-4 ρ ((Derivative[0,0,1][S][r,t,ρ])^2)+8 S[r,t,ρ] ((1-ρ Derivative[0,0,1][B][r,t,ρ]) Derivative[0,0,1][S][r,t,ρ]+ρ Derivative[0,0,2][S][r,t,ρ])-(E^(B[r,t,ρ])) ρ ((S[r,t,ρ])^4) (8+Derivative[0,1,0][B][r,t,ρ] (2 Derivative[1,0,0][B][r,t,ρ]+Derivative[1,0,0][C][r,t,ρ])+Derivative[0,1,0][C][r,t,ρ] (Derivative[1,0,0][B][r,t,ρ]+2 Derivative[1,0,0][C][r,t,ρ]))+((S[r,t,ρ])^2) (2 ρ ((Derivative[0,0,1][B][r,t,ρ])^2)+2 Derivative[0,0,1][C][r,t,ρ]+Derivative[0,0,1][B][r,t,ρ] (-4+ρ Derivative[0,0,1][C][r,t,ρ])+ρ (((Derivative[0,0,1][C][r,t,ρ])^2)-2 Derivative[0,0,2][B][r,t,ρ]-((Derivative[1,0,0][F][r,t,ρ])^2)+24 (E^(B[r,t,ρ])) Derivative[0,1,0][S][r,t,ρ] Derivative[1,0,0][S][r,t,ρ]))),Derivative[1,1,0][F][r,t,ρ]->(1/(2 ρ ((S[r,t,ρ])^2)))(8 ρ Derivative[0,0,1][S][r,t,ρ] Derivative[0,1,0][S][r,t,ρ]+2 ρ S[r,t,ρ] (3 Derivative[0,0,1][S][r,t,ρ] Derivative[0,1,0][B][r,t,ρ]-4 Derivative[0,1,1][S][r,t,ρ]-3 Derivative[0,1,0][S][r,t,ρ] Derivative[1,0,0][F][r,t,ρ]+(Derivative[0,0,1][A][r,t,ρ]+2 Derivative[0,1,0][F][r,t,ρ]) Derivative[1,0,0][S][r,t,ρ])-((S[r,t,ρ])^2) ((-2+2 ρ Derivative[0,0,1][B][r,t,ρ]+ρ Derivative[0,0,1][C][r,t,ρ]) Derivative[0,1,0][B][r,t,ρ]+(2+ρ Derivative[0,0,1][B][r,t,ρ]+2 ρ Derivative[0,0,1][C][r,t,ρ]) Derivative[0,1,0][C][r,t,ρ]+ρ (-2 Derivative[0,1,1][B][r,t,ρ]-(Derivative[0,0,1][A][r,t,ρ]+2 Derivative[0,1,0][F][r,t,ρ]) Derivative[1,0,0][B][r,t,ρ]+Derivative[1,0,0][A][r,t,ρ] Derivative[1,0,0][F][r,t,ρ]+2 Derivative[1,0,1][A][r,t,ρ]))),Derivative[0,2,0][S][r,t,ρ]->(1/(6 ρ ((S[r,t,ρ])^2)))(E^(-B[r,t,ρ])) (-(E^(B[r,t,ρ])) ρ ((S[r,t,ρ])^3) (((Derivative[0,1,0][B][r,t,ρ])^2)+Derivative[0,1,0][B][r,t,ρ] Derivative[0,1,0][C][r,t,ρ]+((Derivative[0,1,0][C][r,t,ρ])^2))+ρ Derivative[0,0,1][S][r,t,ρ] (Derivative[0,0,1][A][r,t,ρ]+2 Derivative[0,1,0][F][r,t,ρ])+S[r,t,ρ] (Derivative[0,0,1][A][r,t,ρ] (1-ρ Derivative[0,0,1][B][r,t,ρ])+ρ Derivative[0,0,2][A][r,t,ρ]-2 (-1+ρ Derivative[0,0,1][B][r,t,ρ]) Derivative[0,1,0][F][r,t,ρ]+2 ρ Derivative[0,1,1][F][r,t,ρ])+3 (E^(B[r,t,ρ])) ρ ((S[r,t,ρ])^2) Derivative[0,1,0][S][r,t,ρ] Derivative[1,0,0][A][r,t,ρ])};

Unprotect[C];

multp={(ρ^(-1/3))(r^(-4)),(ρ^(-2/3)) (r^(-1))((S[r,t,ρ])^2),(ρ^(-4/3))r 12 ((S[r,t,ρ])^3),(ρ^(-2/3))r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(ρ^(-2/3))r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(ρ^(-2/3))r 2 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(r^2)6 ((S[r,t,ρ])^3),(r^2)8 ((S[r,t,ρ])^2)(E^(B[r,t,ρ]))(ρ^(-1/3))};
multp={(ρ^(-1/3))(r^(-4)),(ρ^(-5/3)) (r^(-1))((S[r,t,ρ])^2),(ρ^(-4/3))r 12 ((S[r,t,ρ])^3),(ρ^(-2/3))r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(ρ^(-2/3))r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(ρ^(-2/3)) 2 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(r^2)6 ((S[r,t,ρ])^3),(r^2)8 ((S[r,t,ρ])^2)(E^(B[r,t,ρ]))(ρ^(-1/3))};

multp={(ρ^(-1/3))(r^(-5)),(ρ^(-5/3)) (r^(-2))((S[r,t,ρ])^2),(ρ^(-4/3)) 12 ((S[r,t,ρ])^3),(ρ^(-2/3))r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(ρ^(-2/3))r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ]))/(ρ^2),(ρ^(-2/3)) 2 r ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(r^2)6 ((S[r,t,ρ])^3),(r^2)8 ((S[r,t,ρ])^2)(E^(B[r,t,ρ]))(ρ^(-1/3))}/r;

multp={(r^(-5)),((1+(ρ^2))/(ρ^2))(ρ^1) (r^(-2))((S[r,t,ρ])^2),(ρ^0) 12 ((S[r,t,ρ])^3),(ρ^0)r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),((1+(ρ^2))/(ρ^2))(ρ^2)r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ]))/(ρ^2),(ρ^0) 2 r ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(r^1)6 ((S[r,t,ρ])^3)/ρ,(r^2)8 ((S[r,t,ρ])^2)(E^(B[r,t,ρ]))(ρ^0)}/r;

multp={(r^(-5)),(1/(ρ^2))(ρ^1) (r^(-2))((S[r,t,ρ])^2),(ρ^0) 12 ((S[r,t,ρ])^3),(ρ^0)r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(ρ^2)r 6 ((S[r,t,ρ])^4)(E^(B[r,t,ρ]))/(ρ^2),(ρ^0) 2 r ((S[r,t,ρ])^4)(E^(B[r,t,ρ])),(r^1)6 ((S[r,t,ρ])^3)/ρ,(r^2)8 ((S[r,t,ρ])^2)(E^(B[r,t,ρ]))(ρ^0)}/r;

eqfrinf=eqf/.replnum/.replcoef/.ρ[[i]]->ρ;

If[compute,
EEnew=FullSimplify[multp((#[[1]]-(#[[1]]/.EE)//Expand)&/@EE)];
First[Timing[(eqf=(EEnew/.repldot/.repldρ/.replrtoz//ReleaseHold)/.replfunc/.(((t ρ))^(pow_))->(t^pow)(ρ^pow))]];
ttot=First[Timing[eqfnew=Simplify2/@eqf]];
eqf=eqfnew;
eqf>>"eqfradial3"<>session<>".dat";
ttot,
eqf=<<("eqfradial3"<>session<>".dat");]

(*234.84400000000002`*)

Rest

eqf//LeafCount

Series[#,{z,0,0}]&/@eqf[[1;;-1]]//Simplify//TableForm

(*Series[#,{ρ,0,0}]&/@eqf[[1;;-1]]//Simplify//TableForm*)

solvefor={Sf,Ff,Sdotf,Cdotf,Bdotf,Af};
solveorder={2,2,1,1,1,2};

Timing[Monitor[Do[
(Subscript[As, ii])=D[eqf[[ii]],Derivative[2,0,0][solvefor[[ii]]][z,t,ρ] ]//Simplify2;
(Subscript[Asb, ii])=SeriesCoefficient[Series[(Subscript[As, ii]),{z,0,0}],0]//Simplify2;
(Subscript[Bs, ii])=D[eqf[[ii]],Derivative[1,0,0][solvefor[[ii]]][z,t,ρ]]//Simplify2;
(Subscript[Bsb, ii])=SeriesCoefficient[Series[(Subscript[Bs, ii]),{z,0,0}],0]//Simplify2;
(Subscript[Cs, ii])=D[eqf[[ii]],solvefor[[ii]][z,t,ρ]]//Simplify2;
(Subscript[Csb, ii])=SeriesCoefficient[Series[(Subscript[Cs, ii]),{z,0,0}],0]//Simplify2;
(Subscript[Ss, ii])=eqf[[ii]]/.solvefor[[ii]]->(0&);
(Subscript[Ssb, ii])=SeriesCoefficient[Series[(Subscript[Ss, ii]),{z,0,0}],0]//Simplify2;
If[compute,
(Subscript[Ssbρi, ii])=SeriesCoefficient[Series[(Subscript[Ss, ii])/.replnum/.replcoef/.nums[[nF,a_,b_]]->0/.ksidρ[[i]]->0/.nums[[a_,dρ,b_]]->0/.ρ[[i]]->ρ/.i->1,{ρ,∞,0}],0];
(Subscript[Ssbρ, ii])=SeriesCoefficient[Series[(Subscript[Ss, ii]),{ρ,0,0}],0];];
,{ii,Length[solvefor]}],ii]]

If[compute,
(Subscript[Ssbρ, 5])=(Subscript[Ssbρ, 5])/.Derivative[a_,0,1][Cdotf_][z,t,0]->0/.Derivative[a_,0,3][Bf_][z,t,0]->0/.Derivative[0,1][ξ][t,0]->0/.Derivative[0,3][ξ][t,0]->0//Simplify;
{(Subscript[Ssbρ, 1]),(Subscript[Ssbρ, 2]),(Subscript[Ssbρ, 3]),(Subscript[Ssbρ, 4]),(Subscript[Ssbρ, 5]),(Subscript[Ssbρ, 6])}>>"sourcesr3"<>session<>".dat";
{(Subscript[Ssbρi, 1]),(Subscript[Ssbρi, 2]),(Subscript[Ssbρi, 3]),(Subscript[Ssbρi, 4]),(Subscript[Ssbρi, 5]),(Subscript[Ssbρi, 6])}>>"sourcesri3"<>session<>".dat";,
{(Subscript[Ssbρi, 1]),(Subscript[Ssbρi, 2]),(Subscript[Ssbρi, 3]),(Subscript[Ssbρi, 4]),(Subscript[Ssbρi, 5]),(Subscript[Ssbρi, 6])}=<<("sourcesri3"<>session<>".dat");
{(Subscript[Ssbρ, 1]),(Subscript[Ssbρ, 2]),(Subscript[Ssbρ, 3]),(Subscript[Ssbρ, 4]),(Subscript[Ssbρ, 5]),(Subscript[Ssbρ, 6])}=<<("sourcesr3"<>session<>".dat");];

{(Subscript[Ssbρ, 1]),(Subscript[Ssbρ, 2]),(Subscript[Ssbρ, 3]),(Subscript[Ssbρ, 4]),(Subscript[Ssbρ, 5]),(Subscript[Ssbρ, 6])}//LeafCount

Boundary condition for Sdot &F:

repldt={Derivative[ddz_,1,ddρ_][Ff][z,t,ρ]-> Derivative[ddz,0,ddρ][Fdtf][z,t,ρ],Derivative[ddz_,1,ddρ_][Sdotf][z,t,ρ]-> Derivative[ddz,0,ddρ][Sdotdtf][z,t,ρ]};

SeriesCoefficient[Series[Sdotseries[z,t,ρ],{z,0,2}],2]

SeriesCoefficient[Series[Fseries[z,t,ρ],{z,0,2}],2]/ρ/.replf4/.replcoef

F[z,t,ρ]/.replfunc

updateN:=(
Do[(Subscript[#, ii+10])=(Subscript[#, ii])/.replnumb/.replnumbρ/.replnum/.replcoef;,{ii,Length[solvefor]}]&/@{As,Asb,Bs,Bsb,Cs,Csb,Ss,Ssb,Ssbρ};
If[powers[[nF]]==1,(Subscript[Bsb, 2+10])=0;(Subscript[Csb, 2+10])=1;,(Subscript[Bsb, 2+10])=1;(Subscript[Csb, 2+10])=0;];
(Subscript[Ssb, 2+10])=If[i==1,0,-(1+(ρ^2))SeriesCoefficient[Series[Fseries[z,t,ρ],{z,0,2}],2]/ρ/.replf4/.replcoef//Evaluate];
If[powers[[nSdot]]==1,(Subscript[Bsb, 3+10])=0;(Subscript[Csb, 3+10])=1;,(Subscript[Bsb, 3+10])=1;(Subscript[Csb, 3+10])=0;];
(Subscript[Ssb, 3+10])=Simplify[-SeriesCoefficient[Series[Sdotseries[z,t,ρ],{z,0,2}],2]/.replcoef,Assumptions->{ρ[[i]]>0}];
);
updateN;
constraint=eqf[[8]]/.repldt/.replnum/.replcoef;
constraintb=SeriesCoefficient[Series[eqf[[8]]/.repldt,{z,0,0}],0]/.replnumb/.replnum/.replcoef;
constraintbρ=SeriesCoefficient[Series[eqf[[8]]/.repldt,{ρ,0,0}],0]/.replnumbρ/.replnumb/.replnum/.replcoef;

D[B[z,t,ρ],t]==Bdot[z,t,ρ]-(1/2)A[z,t,ρ]Dr[B,1];
Bfdt=Derivative[0,1,0][Bf][z,t,ρ]/.Solve[D[B[z,t,ρ],t]==Bdot[z,t,ρ]-(1/2)A[z,t,ρ]Dr[B,1]/.replfunc,Derivative[0,1,0][Bf][z,t,ρ]][[1]]//Expand//FullSimplify
Bfdtbr=SeriesCoefficient[Series[Bfdt,{z,0,0}],0]/.replnumb/.replcoef;
Bfdtr=Bfdt/.replnum/.replcoef;
Bfdti=(1/(6 (t^4) z)) (-4-2 t z ((t^2)-t z+(z^2)) Af[z,t,ρ]+6 (t^4) Bdotf[z,t,ρ]+12 (t^4) (z^4) Af[z,t,ρ] Bf[z,t,ρ]+6 (t^4) Cdotf[z,t,ρ]-2 (t^2) (z^2) Af[z,t,ρ] ξ[t,ρ] (-2 t+3 z+3 t z ξ[t,ρ])-2 t ξ[t,ρ] (6+ξ[t,ρ] (7 t+z+t ξ[t,ρ] (4 t+3 z+3 t z ξ[t,ρ])))+4 t z (1+3 t ξ[t,ρ] (1+t ξ[t,ρ])) Derivative[1,0][ξ][t,ρ]+12 (t^4) Bf[z,t,ρ] ((((1+z ξ[t,ρ]))^2)-2 (z^2) Derivative[1,0][ξ][t,ρ])+3 (t^4) (z^5) Af[z,t,ρ] Derivative[1,0,0][Bf][z,t,ρ]+3 (t^4) z ((((1+z ξ[t,ρ]))^2)-2 (z^2) Derivative[1,0][ξ][t,ρ]) Derivative[1,0,0][Bf][z,t,ρ]);
Bfdtir=Bfdti/.replnum/.replcoef/.i->1;
Bfdtzir=SeriesCoefficient[Series[Bfdti,{z,0,0}],0]/.replnumb/.replcoef/.i->1;
Cfdt=Derivative[0,1,0][Cf][z,t,ρ]/.Solve[D[C[z,t,ρ],t]==Cdot[z,t,ρ]-(1/2)A[z,t,ρ]Dr[C,1]/.replfunc,Derivative[0,1,0][Cf][z,t,ρ]][[1]]//Expand//FullSimplify
Cfdtbr=SeriesCoefficient[Series[Cfdt,{z,0,0}],0]/.replnumb/.replcoef;
Cfdtr=Cfdt/.replnum/.replcoef;

a4pρ=SeriesCoefficient[Series[a4p[t,ρ],{ρ,∞,0}],0]/.Derivative[0,1][(Subscript[fs, 4])][t,ρ]->0//Simplify;
a4pr=a4p[t,ρ]/.replcoef
a4pir=a4pρ/.replcoef/.i->1
f4pρ=SeriesCoefficient[Series[f4sp[t,ρ],{ρ,0,0}],0]//Simplify;
f4pr=f4sp[t,ρ]/.replcoef
f4pρr=f4pρ/.replnumbρ/.Derivative[0,2][(Subscript[a, 4])][t,0]->D2Specρ[[nρ]].a4/.Derivative[0,2][(Subscript[b, 4])][t,0]->D2Specρ[[nρ]].b4/.Derivative[0,2][(Subscript[c, 4])][t,0]->D2Specρ[[nρ]].c4

Series[f4sp[tp,ρp],{ρp,0,0}]

Off[Power::infy]
Off[Infinity::indet]
Off[InterpolatingFunction::dmval];
Off[InterpolatingFunction::dmvali];

Evolution procedures

Compute the location of the apparent horizon

ksieq=3 ((S[z,t,ρ])^2) Derivative[0,1,0][S][z,t,ρ]-D[(F[z,t,ρ]+(addξ)'[ρ])S[z,t,ρ](E^(- B[z,t,ρ]))/.z->(1/(1+addξ[ρ])),ρ]+(3/2)(((F[z,t,ρ]+(addξ)'[ρ]))^2)Dr[S,1](E^(- B[z,t,ρ]))/.addξ[ρ]->(1/z)-1;
ksieq=%/.S->(S[#1,#2,#3]((#3)^(1/3))&)/.B->(B[#1,#2,#3]-(2/3)Log[#3]&)/.C->(C[#1,#2,#3]+(4/3)Log[#3]&)//Simplify;
θeq=(ksieq/ρ)/.addξ->(1/AH[#]-1&)//Simplify2
θeqf=Simplify2[θeq/.repldot/.replfunc]//Simplify2;
θeqr=θeqf/.replnum/.replcoef/.replnumtoint/.z->zpos/.zpos->AH[ρ[[i]]]//Simplify;
θeqr=Simplify[θeqr]/.AH[ρ[[i]]]->trial[[i]]/.AH'[ρ[[i]]]->trialdρ[[i]]/.AH''[ρ[[i]]]->trialdρ2[[i]];

θeqinf=SeriesCoefficient[Series[θeqr/.ρ[[i]]->ρ,{ρ,∞,0}],0];

θeqρser=Series[θeqf,{ρ,0,0}]//Simplify;
θeqbρ0=SeriesCoefficient[θeqρser,0];
θeqbρ1=SeriesCoefficient[θeqρser,0];
θeqbρ0r=θeqbρ0/.replnumbρ/.replcoef/.replnumtoint/.z->zpos/.zpos->AH[0]/.AH[0]->trial[[i]]/.AH'[0]->trialdρ[[i]]/.AH''[0]->trialdρ2[[i]]/.i->nρ//Simplify;
θeqbρ1r=θeqbρ1/.replnumbρ/.replcoef/.replnumtoint/.z->zpos/.zpos->AH[0]/.AH[0]->trial[[i]]/.AH'[0]->trialdρ[[i]]/.AH''[0]->trialdρ2[[i]]/.i->nρ//Simplify;

jaceq=D[θeq/.z->zpos/.zpos->AH[ρ],AH[ρ]]//FullSimplify
jaceqr=jaceq/.AH[ρ]->z/.repldot/.replfunc/.replnum/.replcoef/.replnumtoint/.z->zpos/.zpos->AH[ρ[[i]]]/.AH[ρ[[i]]]->trial[[i]]/.AH'[ρ[[i]]]->trialdρ[[i]]/.AH''[ρ[[i]]]->trialdρ2[[i]];
jaceqrdρ=D[θeqr,trialdρ[[i]]]//Simplify;
jaceqrdρ2=D[θeqr,trialdρ2[[i]]]//Simplify;

If[compute,
First[Timing[jaceqρ0=D[θeqbρ0/.z->zpos/.zpos->AH[0],AH[0]]//Simplify]];
ttot=First[Timing[jaceqρ1=D[θeqbρ1/.z->zpos/.zpos->AH[0],AH[0]]//Simplify]];
jaceqρ1>>"jaceq"<>session<>".dat";
ttot,
jaceqρ1=<<("jaceq"<>session<>".dat");]

22.219000000000005`;

jaceqρr=jaceqρ1/.AH[0]->z/.replnumbρ/.replcoef/.replnumtoint/.z->zpos/.zpos->AH[0]/.AH[0]->trial[[i]]/.AH'[0]->trialdρ[[i]]/.AH''[0]->trialdρ2[[i]]/.i->nρ;
jaceqρrdρ=(*D[θeqbρ1r,trialdρ[[nρ]]]//Simplify*)0;
jaceqρrdρ2=D[θeqbρ1r,trialdρ2[[nρ]]]//Simplify;

If[compute,
ttot=First[Timing[jaceqρi=D[θeqinf,trial[[i]]]//Simplify]];
jaceqρi>>"jaceq2"<>session<>".dat";
ttot,
jaceqρi=<<("jaceq2"<>session<>".dat");]

26.921000000000006`;

θeqinfr=SeriesCoefficient[Series[θeqr/.interpolateCheb[nF,a_,i,trial[[i]]]->0/.ksidρ[[i]]->0/.interpolateCheb[a_,dρ,i,trial[[i]]]->0/.ρ[[i]]->ρ,{ρ,∞,0}],0]//Simplify;
jaceqρir=SeriesCoefficient[Series[jaceqr/.interpolateCheb[nF,a_,i,trial[[i]]]->0/.ksidρ[[i]]->0/.interpolateCheb[a_,dρ,i,trial[[i]]]->0/.ρ[[i]]->ρ,{ρ,∞,0}],0];
jaceqρirdρ=SeriesCoefficient[Series[jaceqrdρ/.interpolateCheb[nF,a_,i,trial[[i]]]->0/.ksidρ[[i]]->0/.interpolateCheb[a_,dρ,i,trial[[i]]]->0/.ρ[[i]]->ρ,{ρ,∞,0}],0];
jaceqρirdρ2=SeriesCoefficient[Series[jaceqrdρ2/.interpolateCheb[nF,a_,i,trial[[i]]]->0/.ksidρ[[i]]->0/.interpolateCheb[a_,dρ,i,trial[[i]]]->0/.ρ[[i]]->ρ,{ρ,∞,0}],0];

(*θnew=θnum+jaceq.δtrial+jaceqdy.δtrialdy+jaceqdy2.δtrialdy2=0*)

Cases[#,interpolateCheb[nfunc_,der_,i,trial[[i]]]->{nfunc,der},30]&/@{θeqr,jaceqr,jaceqrdρ,jaceqrdρ2,jaceqρrdρ2};
listused=Flatten[%,1]//Union
%//Length

replintCheb=interpolateCheb[nfunc_,der_,i_,trial[[i_]]]:>(tabCheb[[Position[listused,{nfunc,der}][[1,1]],i]]);
{θeqrtab,jaceqrtab,jaceqrdρtab,jaceqrdρ2tab,jaceqrtabρ,jaceqrdρtabρ,jaceqρrdρ2tab,jaceqρirtab,jaceqρirdρtab,jaceqρirdρ2tab}={θeqr,jaceqr,jaceqrdρ,jaceqrdρ2,jaceqρr,jaceqρrdρ,jaceqρrdρ2,jaceqρir,jaceqρirdρ,jaceqρirdρ2}/.replintCheb;

(*θnew=θnum+jaceq.δtrial+jaceqdy.δtrialdy+jaceqdy2.δtrialdy2=0*)

jacobian4:=(
jnum={jaceqρirtab/.i->1}~Join~Table[jaceqrtab,{i,2,nρ-1}]~Join~{jaceqrtabρ};
jnumdρ={jaceqρirdρtab/.i->1}~Join~Table[jaceqrdρtab,{i,2,nρ-1}]~Join~{jaceqrdρtabρ};
jnumdρ2={jaceqρirdρ2tab/.i->1}~Join~Table[jaceqrdρ2tab,{i,2,nρ-1}]~Join~{jaceqρrdρ2tab};
DiagonalMatrix[jnum]+DiagonalMatrix[jnumdρ].DSpecρ+DiagonalMatrix[jnumdρ2].D2Specρ(*+Table[If[i==nρ+1,DSpecρ[[i,k]],0],{i,nρ},{k,nρ}]*)
)

updateTrial:=(
tabCheb=Table[interpolateCheb[listused[[k,1]],listused[[k,2]],i,trial[[i]]],{k,Length[listused]},{i,nρ}];
trialdρ=DSpecρ.trial;
trialdρ2=D2Specρ.trial;)
computeθ[try_]:=(
trial=try;
updateTrial;
θnum={θeqinfr/.i->1}~Join~Table[θeqrtab,{i,2,nρ-1}]~Join~{ θeqbρ1r})
ϵnew=(10^(-12));
findAHnewton[guess_:1,dostep_:True]:=(
tries={If[guess===1,Table[1,{nρ}],guess]};
If[dostep,step[False,-1,3]];
θ=1;
Monitor[
While[θ>ϵnew,
θ=Max[Abs[computeθ[tries[[-1]]]]];
lp={tries[[-1]]//pl2,θnum//pl2};
If[θ>ϵnew,
AppendTo[tries,trial+LinearSolve[jacobian4,-θnum]]];
],{θ,lp}];
ahpos=tries[[-1]])

Steps

Interpolation[{{0,ksilist[[-5]]},{1,ksilist[[-4]]},{2,ksilist[[-3]]},{3,ksilist[[-2]]},{4,ksilist[[-1]]}},InterpolationOrder->4][5]//Simplify
Derivative[1][Interpolation[{{0,ksilist[[-5]]},{1,ksilist[[-4]]},{2,ksilist[[-3]]},{3,ksilist[[-2]]},{4,ksilist[[-1]]},{5,ksi}},InterpolationOrder->4]][5]//FullSimplify

Integrate[Interpolation[{{0,Bdtlist[[-5]]},{1,Bdtlist[[-4]]},{2,Bdtlist[[-3]]},{3,Bdtlist[[-2]]},{4,Bdtlist[[-1]]}},InterpolationOrder->3][T],{T,4,5}]//Simplify

Table[(Subscript[Ssb, kk+10]),{kk,6}]

Ssbρizi5=SeriesCoefficient[Series[(Subscript[Ssb, 5+10]),{ρ[[i]],∞,0}],0]

step[change_,append_:-1,nsolve_:6,compAH_:False,compConstr_:False]:=(
Do[
nums[[kk+2,1]]=Table[
If[kk<3||kk==6,Anum=SparseArray[(Subscript[As, kk+10])[[1;;-2]]~Join~{(Subscript[Asb, kk+10])}]];
Bnum=SparseArray[(Subscript[Bs, kk+10])[[1;;-2]]~Join~{(Subscript[Bsb, kk+10])}];
Cnum=SparseArray[(Subscript[Cs, kk+10])[[1;;-2]]~Join~{(Subscript[Csb, kk+10])}];
If[i==nρ,
Snum=SparseArray[(Subscript[Ssbρ, kk+10])[[1;;-2]]~Join~{If[kk==5,0,(Subscript[Ssb, kk+10])]}],
If[i==1,
Snum=If[(Subscript[Ssbρi, kk])===0,SparseArray[Table[0,{nz-1}]~Join~{If[kk==5,Ssbρizi5,(Subscript[Ssb, kk+10])]}],SparseArray[(Subscript[Ssbρi, kk])[[1;;-2]]~Join~{If[kk==5,Ssbρizi5,(Subscript[Ssb, kk+10])]}]],
Snum=SparseArray[(Subscript[Ss, kk+10])[[1;;-2]]~Join~{(Subscript[Ssb, kk+10])}]]];
If[kk<3||kk==6,
eqm=DiagonalMatrix[Anum].D2Spec+DiagonalMatrix[Bnum].DSpec+DiagonalMatrix[Cnum],
eqm=DiagonalMatrix[Bnum].DSpec+DiagonalMatrix[Cnum]];
If[i==15,(Subscript[eqmt, kk])=eqm];(* to test invertability*)
LinearSolve[eqm,-Snum]
,{i,nρ}];
(*If[kk+2≠nF,
AppendTo[nums[[kk+2,1]],Table[-Sum[DSpecρ[[nρ,iρ]]nums[[kk+2,1,iρ,iz]],{iρ,nρ-1}]/DSpecρ[[nρ,nρ]],{iz,nz}]],
AppendTo[nums[[kk+2,1]],Table[0,{nz}]]];*)
update[kk+2];
,{kk,1,nsolve}];
If[nsolve==6,
Cdt=Table[filterCheb[Cfdtr[[1;;-2]]~Join~{Cfdtbr}],{i,nρ}];
Cdt=Transpose[filterChebρ[filterBoundDer[#]]&/@Transpose[Cdt]];
(*Bdt={Bfdtir[[1;;-2]]~Join~{Bfdtzir}}~Join~Table[filterCheb[Bfdtr[[1;;-2]]~Join~{Bfdtbr}],{i,2,nρ}];*)
Bdt={Cdt[[1]]}~Join~Table[filterCheb[Bfdtr[[1;;-2]]~Join~{Bfdtbr}],{i,2,nρ-1}]~Join~{Cdt[[-1]]};
Bdt=Transpose[filterChebρ[filterBoundDer[#]]&/@Transpose[Bdt]];

updateCoef;
a4dt=filterChebρ[{a4pir}~Join~Table[a4pr,{i,2,nρ}]];
f4dt=filterChebρ[{0}~Join~Table[If[i==nρ,f4pρr,f4pr],{i,2,nρ}]];];
If[compConstr,
If[Length[Sdotlist]>4&&(timelist[[-5]]-timelist[[-4]])==(timelist[[-2]]-timelist[[-1]]),
nums[[nSdotdt,1]]=((3 Sdotlist[[-4]]-16 Sdotlist[[-3]]+36 Sdotlist[[-2]]-48 Sdotlist[[-1]]+25 nums[[nSdot,1]])/(12 dt));
nums[[nFdt,1]]=((3 Flist[[-4]]-16 Flist[[-3]]+36 Flist[[-2]]-48 Flist[[-1]]+25 nums[[nF,1]])/(12 dt)),
nums[[nSdotdt,1]]=Table[Derivative[1][Interpolation[Table[{timelist[[j]],(Sdotlist~Join~{nums[[nSdot,1]]})[[j,i,k]]},{j,Length[Sdotlist]+1}],InterpolationOrder->Min[order-1,aa-1]]][time],{i,nρ},{k, nz}];
nums[[nFdt,1]]=Table[Derivative[1][Interpolation[Table[{timelist[[j]],(Flist~Join~{nums[[nF,1]]})[[j,i,k]]},{j,Length[Flist]+1}],InterpolationOrder->Min[order-1,aa-1]]][time],{i,nρ},{k, nz}];
];
update[{nSdotdt,nFdt}];
connum=Table[
If[i==nρ,constraintbρ[[1;;-2]]~Join~{constraintb},constraint[[1;;-2]]~Join~{constraintb}],
{i,nρ}];
mcon=N3[connum[[2;;-1,1;;-2]]];
AppendTo[mcons,{time,mcon}]];

If[append≠ -1,
Bdtlist=Append[Drop[Bdtlist,append],Bdt];
Cdtlist=Append[Drop[Cdtlist,append],Cdt];
Sdotlist=Append[Drop[Sdotlist,append],nums[[nSdot,1]]];
Flist=Append[Drop[Flist,append],nums[[nF,1]]];

a4dtlist=Append[Drop[a4dtlist,append],a4dt];
f4dtlist=Append[Drop[f4dtlist,append],f4dt];
];

If[compAH=!=False,ahpos=findAHnewton[compAH]];
If[change,
If[order==4,
nums[[nB,1]]=nums[[nB,1]]+(dt/24)(55 Bdtlist[[-1]]-59 Bdtlist[[-2]]+37 Bdtlist[[-3]]-9 Bdtlist[[-4]]);
nums[[nC,1]]=nums[[nC,1]]+(dt/24)(55 Cdtlist[[-1]]-59 Cdtlist[[-2]]+37 Cdtlist[[-3]]-9 Cdtlist[[-4]]);
a4=a4+(dt/24)(55 a4dtlist[[-1]]-59a4dtlist[[-2]]+37 a4dtlist[[-3]]-9 a4dtlist[[-4]]);
fs4=fs4+(dt/24)(55 f4dtlist[[-1]]-59f4dtlist[[-2]]+37 f4dtlist[[-3]]-9 f4dtlist[[-4]])
,
nums[[nB,1]]=nums[[nB,1]]+(dt/720)(1901 Bdtlist[[-1]]-2774 Bdtlist[[-2]]+2616Bdtlist[[-3]]-1274 Bdtlist[[-4]]+251 Bdtlist[[-5]]);
nums[[nC,1]]=nums[[nC,1]]+(dt/720)(1901 Cdtlist[[-1]]-2774 Cdtlist[[-2]]+2616Cdtlist[[-3]]-1274 Cdtlist[[-4]]+251 Cdtlist[[-5]]);
a4=a4+(dt/720)(1901 a4dtlist[[-1]]-2774 a4dtlist[[-2]]+2616 a4dtlist[[-3]]-1274 a4dtlist[[-4]]+251 a4dtlist[[-5]]);
fs4=fs4+(dt/720)(1901f4dtlist[[-1]]-2774 f4dtlist[[-2]]+2616 f4dtlist[[-3]]-1274f4dtlist[[-4]]+251 f4dtlist[[-5]]);
];
update[{nB,nC}];

time=time+dt;
AppendTo[timelist,time];
updateCoef;
];
)//Timing

dts=Table[(10^(i/2)),{i,-8,-1}]~Join~{1-Sum[(10^(i/2)),{i,-8,-1}]}~Join~Table[1,{points-1}]
startcon=False;
startup[computeAH_]:=Module[{},
timelist={t};
Bdtlist={};
Cdtlist={};
Flist={};
Sdotlist={};
a4dtlist={};
f4dtlist={};
cons={};
mcons={};
Catch[
For[aa=1,aa≤Length[dts],aa++,
step[False,0,6,If[aa==1&&computeAH,1,False],startcon||True];
time=time+dt dts[[aa]];
nums[[nB,1]]=Table[nums[[nB,1,i,k]]+Integrate[Interpolation[Table[{timelist[[j]],Bdtlist[[j,i,k]]},{j,Length[Bdtlist]}],InterpolationOrder->Min[order-1,aa-1]][T],{T,timelist[[-1]],time}],{i,nρ},{k, nz}];
nums[[nC,1]]=Table[nums[[nC,1,i,k]]+Integrate[Interpolation[Table[{timelist[[j]],Cdtlist[[j,i,k]]},{j,Length[Cdtlist]}],InterpolationOrder->Min[order-1,aa-1]][T],{T,timelist[[-1]],time}],{i,nρ},{k, nz}];

If[aa==1&&computeAH,Print["Deviation of AH from 1: ",ahpos-1//Abs//Max]];
If[startcon,
AppendTo[cons,connum//pl]];

a4=Table[a4[[i]]+Integrate[Interpolation[Table[{timelist[[j]],a4dtlist[[j,i]]},{j,Length[a4dtlist]}],InterpolationOrder->Min[order-1,aa-1]][T],{T,timelist[[-1]],time}],{i,nρ}];
fs4=Table[fs4[[i]]+Integrate[Interpolation[Table[{timelist[[j]],f4dtlist[[j,i]]},{j,Length[f4dtlist]}],InterpolationOrder->Min[order-1,aa-1]][T],{T,timelist[[-1]],time}],{i,nρ}];
If[Max[Abs[nums[[nB,1]]]]>(10^4),Throw[Null]];
AppendTo[timelist,time];
update[{nB,nC}];
updateCoef;
];
timelist=timelist[[1;;-6-10]]~Join~timelist[[-7;;-1]];
{Bdtlist,Cdtlist,Sdotlist,Flist}=(#[[1;;-6-9]]~Join~#[[-points+1;;-1]])&/@{Bdtlist,Cdtlist,Sdotlist,Flist};
]];
Off[Power::infy]

Clear[t]

{Derivative[1][Interpolation[{{0,0},{1 ,ksiA},{2,ksiB},{3,ksiC}},InterpolationOrder->3]][0]/((t-tini)/3),
Derivative[2][Interpolation[{{0,0},{1,ksiA},{2,ksiB},{3,ksiC}},InterpolationOrder->3]][0]/((((t-tini)/3))^2),
Derivative[3][Interpolation[{{0,0},{1,ksiA},{2,ksiB},{3,ksiC}},InterpolationOrder->3]][0]/((((t-tini)/3))^3)}//Simplify

initKsi[tini_,profileB_,profileC_,f4in_,a4in_,a4s_,startksi_:1,dorun_:True,grids_:{1.05`60,1.1`60}]:=(
ksiinnum=If[startksi===1,Table[(((a4s[[1]]))^(1/4))-1,{nρ}],{startksi[[1,2]]}~Join~Interpolation[startksi[[2;;-1]]][ρ[[2;;-1]]]];
(*ksiin=Interpolation[Thread[{y,ksiinnum2}]][#2]&;*)
ksiint=Interpolation[Thread[{ρ,ksiinnum}][[2;;-1]]];
ksiinf=ksiinnum[[1]];
ksiinρ=If[#2==∞,ksiinf,ksiint[#2]]&;

Do[
a4IRnow=a4s[[kk]];
run[0,tini,0,profileB,profileC,f4in,(a4in[#]-a4IRnow&),False,ksiinρ,gridsize->grids,execute->False,print->False];
findAHnewton[(*Table[1,{2nρ}]*)];
a4old=a4;
Print["finished ξ computation with regulator ",a4IRnow//N2," after ",Length[tries]," guesses"];
ksiinnum=-(1-(1/(tries[[-1]])))+ksi//filterChebρ2;

ksiint=Interpolation[Thread[{ρ,ksiinnum}][[2;;-1]]];
ksiinf=ksiinnum[[1]];
ksiinρ=If[#2==∞,ksiinf,ksiint[#2]]&;

,{kk,Length[a4s]}];
If[dorun,
dtn=0.0005`60 ((30^2)/((nρ)^2));
run[dtn,tini,10 dtn+tini,profileB,profileC,f4in,(a4in[#]-a4IRnow&),False,ksiinρ,gridsize->{1.025`60,grids[[2]]},execute->True,changegrid->False,print->False,computeAH->True];
Print[t," ",10dtn+tini];
ksiA=-(1-(1/(findAHnewton[])));
Do[step[True,1,6,If[id==10,1,False],True],{id,10}];
ksiB=-(1-(1/(findAHnewton[])));
Do[step[True,1,6,If[id==10,1,False],True],{id,10}];
ksiC=-(1-(1/(findAHnewton[])));
Print["Finished small run till t=",t," with constraint now ",N3[connum[[2;;-1,1;;-2]]]];
ksidts=filterChebρ2/@{ksiinnum,((18 ksiA-9 ksiB+2 ksiC)/(2 t-2 tini)),-((9 (5 ksiA-4 ksiB+ksiC))/(((t-tini))^2)),((27 (3 ksiA-3 ksiB+ksiC))/(((t-tini))^3))};
ksidtinnum=2filterChebρ2[2 ksiA-ksiB/2]/( t-tini);
ksidt2innum=4filterChebρ2[-2ksiA+ksiB]/(((t-tini))^2),
ksidtinnum=Table[0,{nρ}]];

ksiint=Interpolation[Thread[{ρ,ksiinnum}][[2;;-1]]];
ksiindt=Interpolation[Thread[{ρ,ksidtinnum}][[2;;-1]]];
ksiininf=ksiinnum[[1]]+#1 ksidtinnum[[1]]&;

ksiin=If[#2==∞,ksiininf[#1],ksiint[#2]+#1 ksiindt[#2]]&;
)

cln

2 ksiA-ksiB/2

Derivative[1][Interpolation[{{0,ksiinnum},{1,ksiA},{2,ksiB}},InterpolationOrder->2]][6]/(t-tini)//Simplify
Derivative[2][Interpolation[{{0,ksiinnum},{1,ksiA},{2,ksiB}},InterpolationOrder->2]][6]/(((t-tini))^2)//Simplify

Derivative[1][Interpolation[{{0,ksi1},{1,ksia},{2,ksib}},InterpolationOrder->3]][2]//Simplify
Derivative[2][Interpolation[{{0,ksi1},{1,ksia},{2,ksib}},InterpolationOrder->3]][2]/(tno-tini)//Simplify

cln

Clear[ksilist]

Derivative[1][Interpolation[{{0,ksilist[[-5]]},{1,ksilist[[-4]]},{2,ksilist[[-3]]},{3,ksilist[[-2]]},{4,ksilist[[-1]]}},InterpolationOrder->3]][4]//Simplify
Derivative[2][Interpolation[{{0,ksilist[[-5]]},{1,ksilist[[-4]]},{2,ksilist[[-3]]},{3,ksilist[[-2]]},{4,ksilist[[-1]]}},InterpolationOrder->3]][4]//Simplify

filterChebρ2[(1/(3updateAHtime dt))(40(*20*) ksidif+10 ksidtdif+2 ksidt2dif)]

updateKsi[(Derivative[2][Function[tp,Interpolation[ksinew[[-4;;-1]]][tp]]][t]+8 Derivative[1][Function[tp,Interpolation[ksinew[[-4;;-1]]][tp]]][t]+300 (ksinew[[-1,2]])/2),20 dt];

+(1/6)(targetksidt/difference)(((#-ti0))^3)

20 0.0005

updateKsi:=(
ksidif=ksilist[[-1]];
ksidtdif=(1/(updateAHtime dt)) ((1/6) (-2 ksilist[[-4]]+9 ksilist[[-3]]-18 ksilist[[-2]]+11 ksilist[[-1]]));
ksidt2difold=ksidt2dif;
ksidt2dif=(1/(((updateAHtime dt))^2)) (-ksilist[[-4]]+4 ksilist[[-3]]-5 ksilist[[-2]]+2 ksilist[[-1]]);
ksiatT=ksi;
ksidtatT=ksidt;
ksidt2atT=Derivative[2,0][ksifunc][t,#]&/@ρ;
tnow=t;
(*ksidt3=filterChebρ2[(1/(3updateAHtime dt))(40(*20*) ksidif+10 ksidtdif+2 ksidt2dif)];*)
If[dif=={},
ksidt3=filterChebρ2[(1/(updateAHtime dt))(150(*20*) ksidif+12 ksidtdif+ksidt2dif)],
ksidt3=filterChebρ2[(1/( updateAHtime dt))(150(*20*) ksidif+8 ksidtdif+4ksidt2dif/4)]];
ksidt3=ksidt3/Max[{Max[Abs[ksidt3]]/Max[{10,70-(t-3/10)60 12/7}],1}];
ksiintfunc=Interpolation[Thread[{ρ,ksiatT}][[2;;-1]]][#2]+(#1-tnow)Interpolation[Thread[{ρ,ksidtatT}][[2;;-1]]][#2]+(1/2)(((#1-tnow))^2)Interpolation[Thread[{ρ,ksidt2atT}][[2;;-1]]][#2]+(1/6)(((#1-tnow))^3)Interpolation[Thread[{ρ,ksidt3}][[2;;-1]]][#2]&;
ksiatinf=ksiatT[[1]]+(#1-tnow)ksidtatT[[1]]+(1/2)(((#1-tnow))^2)ksidt2atT[[1]]+(1/6)(((#1-tnow))^3)ksidt3[[1]]&;

ksifunc=If[#2==∞,ksiatinf[#1],ksiintfunc[#1,#2]]&;
ksidtfunc=Derivative[1,0][ksifunc];
ksilist=Drop[ksilist,1];
dif={ksidif//pl2, ksidtdif//pl2,ksidt2dif//pl2,ksidt3//pl2};
);

Helpful quantities

cl

(1/dt)Derivative[1][Interpolation[{{0,Sdotlist[[-5]]},{1,Sdotlist[[-4]]},{2,Sdotlist[[-3]]},{3,Sdotlist[[-2]]},{4,Sdotlist[[-1]]}},InterpolationOrder->4]][5]//Simplify
%/.Sdotlist->Flist;

plotConstraint:=(
step[False,-1,6,False,True];
connum//pl)

(1/dt)Derivative[1][Interpolation[{{0,Sdotlist[[-5]]},{1,Sdotlist[[-4]]},{2,Sdotlist[[-3]]},{3,Sdotlist[[-2]]},{4,Sdotlist[[-1]]}},InterpolationOrder->4]][4]//Simplify
%/.Sdotlist->Flist

savestate:={t,ii,nums[[nB,1]],a4,f4,ksi,ksidt,ksilist,dt,tend,fixAH}

loadstate[file_]:=(
{t,ii,nums[[nB,1]],a4,f4,ksi,ksidt,ksilist,dt,tend,fixAH}=<<file;
run;)

eqsinB=0==eqf[[1]]/.ξ->(0&)/.Sf->(Sf[#1]&)/.Bf->(Bf[#1]&)/.z->zp//Simplify

eventsinB=S[z,t,y]==0/.replfunc/.ξ->(0&)/.Sf->(Sf[#1]&)/.z->zp

singularityfromBaty[Bf_]:=(
ndsol=NDSolve[{16 ((Bf[zp])^2) ((zp^3)+(zp^8) Sf[zp])+8 (zp^4) Bf[zp] (1+(zp^5) Sf[zp]) Bf'[zp]+Sf[zp] (40+(zp^10) ((Bf'[zp])^2))+zp ((zp^4) ((Bf'[zp])^2)+20 Sf'[zp]+2 zp Sf''[zp])==0,Sf[0.005]==0,Sf'[0.005]==0},Sf,{zp,0.005,8.},Method->{"EventLocator", "Event":>(1/zp)+(zp^4) Sf[zp]>0},AccuracyGoal->15,PrecisionGoal->15][[1]];
low=InterpolatingFunctionDomain[First[Sf/.ndsol[[1]]]][[1,-1,-1]];
If[machinep,{low,(1/low)+(low^4) Sf[low]/.ndsol},SetPrecision[{low,Sf[low]/.ndsol},30]]);

singularityfromBaty[B_]:=(
ndsol=NDSolve[{s''[zp]+2 s'[zp]/zp+(1/2)((B'[zp])^2) s[zp]==0,s[0.001]==1000,s'[0.001]==-1/(0.001^2)},s,{zp,0.001,8.},Method->{"EventLocator", "Event":>s[zp]>0}][[1]];
low=InterpolatingFunctionDomain[First[s/.ndsol[[1]]]][[1,-1,-1]];
If[machinep,{low,s[low]/.ndsol},SetPrecision[{low,s[low]/.ndsol},30]]);

singularityfromB[B_]:=(
singularityfromBaty[Function[zp,B[zp,#]]]&/@y)

S[z,t,ρ]/.replfunc;
Sr=%/.replnum/.replcoef;
Sri=Sr/.replnumtoint/.z->zpos/.zpos->ahpos[[i]];
AHarea:=Table[(Sri^3),{i,nρ}];

Sri//Definition

AHarea//Definition

Run

Options[run]={print->True,gridsize->{1.01`60,4.`60},onlyrun->8,execute->True,changegrid->False,computeAH->False};
run[timestep_,tini_,end_,profileB_,profileC_,f4in_,a4in_,fixah_,ksiin_:(0#&),OptionsPattern[]]:=(
dt=timestep;
tend=end;
fixAH=fixah;
If[fixah&&(ksiin===(0#&)),ksiini=initKsi[profileB,profileC,f4in,a4in,{2,1.5,1,0.5,0.4,0.3}],ksiini=ksiin];
If[OptionValue[onlyrun]==8,
init[tini,profileB,profileC,f4in,a4in,ksiini,OptionValue[gridsize][[1]],OptionValue[gridsize][[2]]];
];
updateAHtime=Round[0.007`60/dt];
If [OptionValue[print],Print["Initialized evolution of B=",profileB["z","ρ"]," and C=",profileC["z","ρ"]," running from ",time," to ",end,". Using (Subscript[n, z]) = ",nz," and (Subscript[n, ρ]) = ",nρ,". Grid setting: ",OptionValue[changegrid]," filter z&ρ: ",nfilterz," ",nfilterρ," ",nfilterρ2," sessionID: ",session," ksiID:",sesksi]];
If[OptionValue[execute],
Catch[
Monitor[
If[OptionValue[onlyrun]==8,
Bfins={Table[{z[[i]],time,ρ[[k]],nums[[nB,1,k,i]]},{k,nρ},{i,nz}]};
Cfins={Table[{z[[i]],time,ρ[[k]],nums[[nC,1,k,i]]},{k,nρ},{i,nz}]};
a4s={Table[{time,ρ[[k]],a4[[k]]},{k,nρ}]};
f4s={Table[{time,ρ[[k]],fs4[[k]]},{k,nρ}]};
b4s={Table[{time,ρ[[k]],b4[[k]]},{k,nρ}]};
c4s={Table[{time,ρ[[k]],c4[[k]]},{k,nρ}]};
ksis={Table[{time,ρ[[k]],ksi[[k]]},{k,nρ}]};
dif={}; 
AHA=AHarea;
AHAs={Table[{time,ρ[[k]],AHA[[k]]},{k,nρ}]};
ttotal=First[Timing[startup[OptionValue[computeAH]]]];
ksilist={};
ahpos={1};
tries={};
mcon=0;];
For[ii=OptionValue[onlyrun],ii≤Round[(end-ti)/timestep],ii++,
If[Max[Abs[nums[[nB,1]]]]>(10^4)||mcon>(10^3),Throw[Null]];
If[Mod[ii-8,25]==0,(*save metric for later use*)
AppendTo[Bfins,Table[{z[[i]],time,ρ[[k]],nums[[nB,1,k,i]]},{k,nρ},{i,nz}]];
AppendTo[Cfins,Table[{z[[i]],time,ρ[[k]],nums[[nC,1,k,i]]},{k,nρ},{i,nz}]];
AppendTo[a4s,Table[{time,ρ[[k]],a4[[k]]},{k,nρ}]];
AppendTo[b4s,Table[{time,ρ[[k]],b4[[k]]},{k,nρ}]];
AppendTo[c4s,Table[{time,ρ[[k]],c4[[k]]},{k,nρ}]];
AppendTo[f4s,Table[{time,ρ[[k]],fs4[[k]]},{k,nρ}]];
AppendTo[ksis,Table[{time,ρ[[k]],ksi[[k]]},{k,nρ}]];
AppendTo[AHAs,Table[{time,ρ[[k]],AHA[[k]]},{k,nρ}]];
(*AHA=Table[(Sri^3),{i,nρ}];*)
];
findAH=If[fixah&&Mod[ii-8,updateAHtime]==0,If[Length[tries]==0,Table[1,{nρ}],tries[[-1]]],False];
printnow=((ii==25)||(ii==15)||(Mod[ii,40]==0))&&OptionValue[print];
If[Mod[ii,1000]==0,{Bfins,Cfins,a4s,b4s,c4s,f4s,ksis,mcons,L,powercoord}>>"result"<>sesksi];

ttotal=ttotal+step[True,1,6,findAH,printnow||True][[1]];

(*AppendTo[cons,pl[connum]];*)

If[findAH=!=False,(*determine AH for fixing ξ*)
AHA=AHarea;
If[OptionValue[changegrid]=="changegrid",
changeGrid[-(1-1/filterChebρ[ahpos])];];
If[OptionValue[changegrid]=="changeksi",
AppendTo[ksilist,-(1-(1/(filterChebρ[ahpos])))];
If[Length[ksilist]==4,updateKsi]];
];
If[printnow,
con=pl[connum];
pla4={pllist[a4,a4in],pllist[b4,b4ini],pllist[c4,c4ini],pllist[fs4,f4in],pllist[ksi,ksiini[ti,#1]&]};
If [OptionValue[print],Print[ii," ",N3[connum[[2;;-1,1;;-2]]]," ",N3[connum[[1;;-1,-1]]]," ",N3[ahpos-1]," ",N3[ksidt3]," ",N3[ksi],"  ",N3[fs4],"  ",N3[b4],"  ",N3[a4],"  ",N2[time]," ",ttotal]];
];
];
,{pla4,lp,con,dif}]],If[OptionValue[computeAH],Print["Deviation of AH from 1: ",findAHnewton[]-1//Abs//Max]]]
);

Plots and helpful quantities

Solve[L2 (((ρCheb2)/(((1-((ρCheb2)^2)))^(1/4))))==x2,ρCheb2]

(Sqrt[-((x2^4)/(2 (L2^4)))+((Sqrt[4 (L2^4) (x2^4)+(x2^8)])/(2 (L2^4)))])/.x2->Slot[2]/.L2->L

intlist[list_]:=Module[{intfunc=Interpolation[Flatten[Table[Thread[{Transpose[list[[i]]][[1]],ρCheb[[1;;nρ]],Transpose[list[[i]]][[3]]}],{i,Length[list]}],1],InterpolationOrder->2,Method-> "Spline"]},
If[#2==∞,intfunc[#1,1],intfunc[#1,(Sqrt[-(((#2)^4)/32)+(1/32) (Sqrt[64 ((#2)^4)+((#2)^8)])])]]&]

intlist[list_]:=Module[{intfunc=Interpolation[Flatten[Table[Thread[{Transpose[list[[i]]][[1]],ρCheb[[1;;nρ]],Transpose[list[[i]]][[3]]}],{i,Length[list]}],1],InterpolationOrder->2,Method-> "Spline"]},
Function[{tp,ρp},If[ρp==∞,intfunc[tp,1],intfunc[tp,(Sqrt[Root[-((ρp)^powercoord)+((ρp)^powercoord) #1+(L^powercoord) ((#1)^(powercoord/2))&,2]])]]]]

intlist[list_]:=Module[{intfunc=Interpolation[Flatten[Table[Thread[{Transpose[list[[i]]][[1]],ρCheb[[1;;nρ]],Transpose[list[[i]]][[3]]}],{i,Length[list]}],1],InterpolationOrder->2,Method-> "Spline"]},
Function[{tp,ρp},If[ρp==∞,intfunc[tp,1],intfunc[tp,coordtrans[ρp]]]]]

Root[-((ρx)^6)+((ρx)^6) #1+(Lx^6) ((#1)^(6/2))&,1]//ToRadicals

(Sqrt[Root[-((ρx)^18)+((ρx)^18) #1+(Lx^18) ((#1)^9)&,1]])/.ρx->.3/.Lx->22

Solve[ρx==Lx (((((ρChebx)/(((1-((ρChebx)^2)))^(1/17)))))^1),ρChebx]

Solve[ρx==Lx (((((ρChebx)/(((1-((ρChebx)^2)))^(1/18)))))^1),ρChebx]

v[t_,ρ_]:=((-1)/(12 (t^4) ρ f4si[t,ρ]))(-1-6 (t^4) a4si[t,ρ]+6 (t^4) b4si[t,ρ]-4 t ksii[t,ρ]-6 (t^2) ((ksii[t,ρ])^2)-4 (t^3) ((ksii[t,ρ])^3)-(Sqrt[-144 (t^8) (ρ^2) ((f4si[t,ρ])^2)+(((1+6 (t^4) a4si[t,ρ]-6 (t^4) b4si[t,ρ]+4 t ksii[t,ρ]+6 (t^2) ((ksii[t,ρ])^2)+4 (t^3) ((ksii[t,ρ])^3)))^2)]))
v2[t_,ρ_]:=((-1)/(12 (t^4) ρ f4si[t,ρ]))(-1-6 (t^4) a4si[t,ρ]+6 (t^4) b4si[t,ρ]-4 t ksii[t,ρ]-6 (t^2) ((ksii[t,ρ])^2)-4 (t^3) ((ksii[t,ρ])^3)+(Sqrt[-144 (t^8) (ρ^2) ((f4si[t,ρ])^2)+(((1+6 (t^4) a4si[t,ρ]-6 (t^4) b4si[t,ρ]+4 t ksii[t,ρ]+6 (t^2) ((ksii[t,ρ])^2)+4 (t^3) ((ksii[t,ρ])^3)))^2)]))
v3[t_,ρ_]:=(If[ρ==∞||ρ==0||t<ti+0.001,0,vtmp1=v[t,ρ];If[Abs[vtmp1]≤1,vtmp1,v2[t,ρ]]])
Elocal[t_,ρ_]:=(1/(12 (t^4)))(1-3 (t^4) a4si[t,ρ]-6 (t^4) b4si[t,ρ]+4 t ksii[t,ρ]+6 (t^2) ((ksii[t,ρ])^2)+4 (t^3) ((ksii[t,ρ])^3)+(Sqrt[-144 (t^8) (ρ^2) ((f4si[t,ρ])^2)+(((1+6 (t^4) a4si[t,ρ]-6 (t^4) b4si[t,ρ]+4 t ksii[t,ρ]+6 (t^2) ((ksii[t,ρ])^2)+4 (t^3) ((ksii[t,ρ])^3)))^2)]));
Plocalρ[t_,ρ_]:=(1/(12 (t^4)))(-1+3 (t^4) a4si[t,ρ]+6 (t^4) b4si[t,ρ]-4 t ksii[t,ρ]-6 (t^2) ((ksii[t,ρ])^2)-4 (t^3) ((ksii[t,ρ])^3)+(Sqrt[-144 (t^8) (ρ^2) ((f4si[t,ρ])^2)+(((1+6 (t^4) a4si[t,ρ]-6 (t^4) b4si[t,ρ]+4 t ksii[t,ρ]+6 (t^2) ((ksii[t,ρ])^2)+4 (t^3) ((ksii[t,ρ])^3)))^2)]));
Plocalθ[t_,ρ_]:=-((2+3 (t^4) a4si[t,ρ]-12 (t^4) c4si[t,ρ]+8 t ksii[t,ρ]+12 (t^2) ((ksii[t,ρ])^2)+8 (t^3) ((ksii[t,ρ])^3))/(12 (t^4)));
Plocaly[t_,ρ_]:=-(1/4) a4si[t,ρ]+((1-3 (t^4) b4si[t,ρ]-3 (t^4) c4si[t,ρ]+4 t ksii[t,ρ]+6 (t^2) ((ksii[t,ρ])^2)+4 (t^3) ((ksii[t,ρ])^3))/(3 (t^4)));

Teffa4[a4_]:=(((-a4))^(1/4))/π;
Teff[e_]:=(((4e/3))^(1/4))/π;

{Thydroρρ[τ_,ρ_],Thydroθθ[τ_,ρ_],Thydroyy[τ_,ρ_]}:={((Elocal[τ,ρ] (1+3 ((vloc[τ,ρ])^2)))/(3-3 ((vloc[τ,ρ])^2)))-((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (ρ ((vloc[τ,ρ])^2)+τ ((vloc[τ,ρ])^3)+ρ (-1+2 τ Derivative[0,1][vloc][τ,ρ])+τ vloc[τ,ρ] (-1+2 ρ Derivative[1,0][vloc][τ,ρ])))/(3 ρ τ (((1-((vloc[τ,ρ])^2)))^(5/2)))),((ρ Elocal[τ,ρ]+((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (ρ-ρ ((vloc[τ,ρ])^2)+2 τ ((vloc[τ,ρ])^3)+ρ τ Derivative[0,1][vloc][τ,ρ]+τ vloc[τ,ρ] (-2+ρ Derivative[1,0][vloc][τ,ρ])))/(τ (((1-((vloc[τ,ρ])^2)))^(3/2)))))/(3 ρ)),((τ Elocal[τ,ρ]+((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (2 ρ ((vloc[τ,ρ])^2)-τ ((vloc[τ,ρ])^3)+ρ (-2+τ Derivative[0,1][vloc][τ,ρ])+vloc[τ,ρ] (τ+ρ τ Derivative[1,0][vloc][τ,ρ])))/(ρ (((1-((vloc[τ,ρ])^2)))^(3/2)))))/(3 τ))};

{Tρρ[τ_,ρ_],Tθθ[τ_,ρ_],Tyy[τ_,ρ_]}:={-(1/(6 (τ^4)))-(1/4) a4si[τ,ρ]+b4si[τ,ρ]-((2 ksii[τ,ρ])/(3 (τ^3)))-(((ksii[τ,ρ])^2)/(τ^2))-((2 ((ksii[τ,ρ])^3))/(3 τ)),-((2+3 (τ^4) a4si[τ,ρ]-12 (τ^4) c4si[τ,ρ]+8 τ ksii[τ,ρ]+12 (τ^2) ((ksii[τ,ρ])^2)+8 (τ^3) ((ksii[τ,ρ])^3))/(12 (τ^4))),-(1/4) a4si[τ,ρ]+((1-3 (τ^4) b4si[τ,ρ]-3 (τ^4) c4si[τ,ρ]+4 τ ksii[τ,ρ]+6 (τ^2) ((ksii[τ,ρ])^2)+4 (τ^3) ((ksii[τ,ρ])^3))/(3 (τ^4)))};

{Thydroρρη[τ_,ρ_],Thydroθθη[τ_,ρ_],Thydroyyη[τ_,ρ_]}:={-((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (ρ ((vloc[τ,ρ])^2)+τ ((vloc[τ,ρ])^3)+ρ (-1+2 τ Derivative[0,1][vloc][τ,ρ])+τ vloc[τ,ρ] (-1+2 ρ Derivative[1,0][vloc][τ,ρ])))/(3 ρ τ (((1-((vloc[τ,ρ])^2)))^(5/2)))),(((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (ρ-ρ ((vloc[τ,ρ])^2)+2 τ ((vloc[τ,ρ])^3)+ρ τ Derivative[0,1][vloc][τ,ρ]+τ vloc[τ,ρ] (-2+ρ Derivative[1,0][vloc][τ,ρ])))/(τ (((1-((vloc[τ,ρ])^2)))^(3/2))))/(3 ρ)),(((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (2 ρ ((vloc[τ,ρ])^2)-τ ((vloc[τ,ρ])^3)+ρ (-2+τ Derivative[0,1][vloc][τ,ρ])+vloc[τ,ρ] (τ+ρ τ Derivative[1,0][vloc][τ,ρ])))/(ρ (((1-((vloc[τ,ρ])^2)))^(3/2))))/(3 τ))};

-((2 (((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))))  (ρ ((vloc[τ,ρ])^2)+τ ((vloc[τ,ρ])^3)+ρ (-1+2 τ Derivative[0,1][vloc][τ,ρ])+τ vloc[τ,ρ] (-1+2 ρ Derivative[1,0][vloc][τ,ρ])))/(3 ρ τ (((1-((vloc[τ,ρ])^2)))^(5/2))))

Thydroρρη//Definition

η:=(((Elocal[τ,ρ])^(3/4))/((Sqrt[2]) (3^(3/4))));

findfactor:=(
etab=Table[{tp,Elocal[tp,0.02]},{tp,ti,t,0.01}];
easym[tp_]=Fit[etab[[-40;;-1]],{(tp^(-4/3)),(tp^(-6/3))},tp];
factor=fac/.FindRoot[Teff[easym[0.6 /fac]]/fac==2.564271967629648`,{fac,0.3}])

findfactor:=(
etab=Table[{tp,((-27)/(5.068 8 (π^2)))a4si[tp,0.02]},{tp,ti,t,0.01}];
easym[tp_]=Fit[etab[[-40;;-1]],{(tp^(-4/3)),(tp^(-6/3))},tp];
factor=fac/.FindRoot[easym[0.6 /fac]/(fac^4)==187,{fac,0.3}])

plotfactor:=Show[Plot[easym[tp],{tp,ti,0.65/factor}],etab//ListPlot,Graphics[Line[{{0.6/factor,0},{0.6/factor,1000}}]]]

findtiso:=(
tthermal=FindRoot[((Tρρ[tp,0.02]-Thydroρρ[tp,0.02])/(Elocal[tp,0.02]))==-0.02,{tp,.4}];
{tp/.(thermal=FindRoot[((Tyy[tp,0.02]-Thydroyy[tp,0.02])/(Elocal[tp,0.02]))==0.02,{tp,.4}]),
ethermal=Elocal[tp,0.02]/.tthermal,
Teff[ethermal],
wiso=tp Teff[ethermal]/.tthermal});

plotv:=Plot3D[v3[tp/factor,yp/factor],{tp,ti factor,t factor},{yp,0.05 factor,factor Min[{ρ[[2]]1.05,(ymax+L/4)}]},ColorFunction->"Rainbow",PlotRange->All,AxesLabel->{Text["τ(fm/c)",BaseStyle->{FontSize->30}],Text["ρ (fm)",BaseStyle->{FontSize->30}],Text["v",BaseStyle->{Italic,FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->25,FaceGrids->{{{0,1,0},{{0,0.2,0.3,0.4,0.5,0.6,0.7},Automatic}},{{-1,0,0},{If[ymax>5,{-0.01,2.5,5,7.5,10},{-0.01,0.25,0.5,0.75,1.0,1.25,1.5}],Automatic}}},FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plotvE:=Plot3D[v3[tp/factor,yp/factor](((-27)/(5.068 8 (π^2)))a4si[tp/factor,yp/ factor](factor^(-4))),{tp,ti factor,t factor},{yp,0.05 factor,factor Min[{ρ[[2]]1.05,(ymax+L/4)}]},ColorFunction->"Rainbow",PlotRange->All,AxesLabel->{Text["τ(fm)",BaseStyle->{FontSize->30}],Text["ρ (fm)",BaseStyle->{FontSize->30}],Text["ε (Subscript[v, ρ]) (GeV/(fm^3))",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->60,FaceGrids->{{{0,1,0},{{0,0.2,0.3,0.4,0.5,0.6,0.7},Automatic}},{{-1,0,0},{If[ymax>5,{-0.01,2.5,5,7.5,10},{-0.01,0.25,0.5,0.75,1.0,1.25,1.5}],Automatic}}},FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plota:=Plot3D[(1/factor)Derivative[1,0][v][tp/factor,yp/factor+(10^(-5))]//Abs,{tp,factor ti+0.01,t factor},{yp,0.0 factor,factor Min[{ρ[[2]]1.05,(ymax+L/4)}]},ColorFunction->"Rainbow",PlotRange->Automatic,AxesLabel->{Text["τ (fm)",BaseStyle->{FontSize->30}],Text["ρ (fm)",BaseStyle->{FontSize->30}],Text["(Subscript[a, ρ]) ((fm^(-1)))",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->50,FaceGrids->{{{0,1,0},{{0,0.2,0.3,0.4,0.5,0.6,0.7},Automatic}},{{-1,0,0},{If[ymax>5,{-0.01,2.5,5,7.5,10},{-0.01,0.25,0.5,0.75,1.0,1.25,1.5}],Automatic}}},FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plotf4:=(mini[tp_]:=yp/.FindMinimum[f4si[tp,yp] yp,{yp,1.2L,0,3L}][[2]];
border[tp_]:=yp/.FindRoot[f4si[tp,yp] yp+0.02,{yp,1.4L,L/2,2L}][[1]];
minis=Table[{tp factor,factor mini[tp]},{tp,0.2+0.01,t+0.01,0.025}];
borders=Table[{tp factor,factor border[tp]},{tp,0.3,t+0.01,0.05}];
plmini=ListLinePlot[{Transpose[Reverse[Transpose[minis]]],Transpose[Reverse[Transpose[borders]]]},PlotStyle->{{Dashed,Black,Thickness[0.01]},{Blue,Thickness[0.007]}}];
Show[ContourPlot[f4si[tp/factor,yp/factor] yp/factor,{yp,0,ymax factor},{tp,ti factor,t factor},ColorFunction->"Rainbow",PlotRange->{{0,ymax factor},{ti factor,t factor},All},FrameLabel->{Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text["τ(fm/c)",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->50,AxesStyle->Thickness[0.004],Mesh->10,Contours->20],plmini])

plotE:=Plot3D[(((-3)/4)(a4si[tp/factor,yp/factor]))tp,{tp,ti factor,t factor},{yp,0,ymax factor},ColorFunction->"Rainbow",PlotRange->All,AxesLabel->{Text["τ(fm/c)",BaseStyle->{FontSize->30}],Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text["ℰ τ",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->50,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plotTlocal:=Plot3D[(1/(factor 5.0677))(((((4(Elocal[tp/factor,yp/factor]-Elocal[tp/factor,ρ[[2]]]+0.002))/(3(π^4)))))^(1/4))(tp^(0/3)),{tp,ti factor,t factor},{yp,0.02,factor(ymax+3)},ColorFunction->"Rainbow",PlotRange->All,AxesLabel->{Text["τ(fm/c)",BaseStyle->{FontSize->30}],Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text["(Subscript[T, eff, local])(GeV)",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->25,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plotT:=Plot3D[(1/(factor 5.0677))(((-(((a4si[tp/factor,yp/factor]-a4si[tp/factor,ρ[[2]]]+0.002))/(π^4))))^(1/4))(tp^(0/3)),{tp,ti factor,t factor},{yp,0.02,factor(ymax+3)},ColorFunction->"Rainbow",PlotRange->All,AxesLabel->{Text["τ(fm/c)",BaseStyle->{FontSize->30}],Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text["(Subscript[T, eff])(GeV)",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->25,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plotPminE:=Plot3D[((Plocalρ[tp/factor,yp/factor]-Elocal[tp/factor,yp/factor]/3)/(Elocal[tp/factor,yp/factor])),{tp,ti factor,t factor},{yp,0.02,factor(ymax+5)},ColorFunction->"Rainbow",PlotRange->All,AxesLabel->{Text["τ(fm/c)",BaseStyle->{FontSize->30}],Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text["(((Subscript[𝒫, ρ]) - ℰ/3)/ℰ)",BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->25,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

plotΔP:=Plot3D[#[[1]][tp/factor,yp/factor],{tp,ti factor,t factor},{yp,0.02,factor(ymax+3)},ColorFunction->"Rainbow",PlotRange->{All,All,{-0.1,0.4}},AxesLabel->{Text["τ(fm)",BaseStyle->{FontSize->30}],Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text[#[[2]],BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->25,FaceGrids->{{{0,1,0},{{0,0.2,0.3,0.4,0.5},Automatic}},{{-1,0,0},{{-0.01,2.5,5,7.5,10},Automatic}}},FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]&/@Transpose[{{((Tρρ[#1,#2]-Thydroρρ[#1,#2])/(Elocal[#1,#2]))&,((Tθθ[#1,#2]-Thydroθθ[#1,#2])/(Elocal[#1,#2]))&,((Tyy[#1,#2]-Thydroyy[#1,#2])/(Elocal[#1,#2]))&},{"(((Subscript[T, ρρ]) - (Subscript[((T^hydro)), ρρ]))/ε)","(((Subscript[T, θθ]) - (Subscript[((T^hydro)), θθ]))/ε)","(((Subscript[T, yy]) - (Subscript[((T^hydro)), yy]))/ε)"}}]

plothydrorho[range_:{-0.9,0.9}]:=Show[Plot3D[-3 #[[1]][tp/factor,yp/factor],{tp,ti factor,t factor},{yp,0.02,factor ρ[[2]]1.05},ColorFunction->Function[{x,y,z},ColorData["Rainbow"][z/Max[Abs[range]]//Abs]],ColorFunctionScaling->False,PlotRange->{All,All,range},AxesLabel->{Text["τ (fm)",BaseStyle->{FontSize->30}],Text["ρ (fm)",BaseStyle->{FontSize->30}],Text[#[[2]],BaseStyle->{FontSize->34}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->50,FaceGrids->{{{0,1,0},{{0,0.2,0.3,0.4,0.5,0.6,0.7},Automatic}},{{-1,0,0},{If[ymax>5,{-0.01,2.5,5,7.5,10},{-0.01,0.25,0.5,0.75,1.0,1.25,1.5}],Automatic}}},FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]&[{((Tρρ[#1,#2]-Thydroρρ[#1,#2])/(Elocal[#1,#2]))&,"(((Subscript[((T^hydro)), ρ]^ρ) - (Subscript[p, ρ]))/(ε/3))"}],
ParametricPlot3D[{{tp,0.001,((3 Thydroρρη[tp/factor,0.001])/(Elocal[tp/factor,0.001]))},{tp,ρ[[2]]1.05 factor/2,((3 Thydroρρη[tp/factor,ρ[[2]]1.05 /2])/(Elocal[tp/factor,ρ[[2]]1.05 /2]))},{tp,ρ[[2]]1.05 factor,((3 Thydroρρη[tp/factor,ρ[[2]]1.05 ])/(Elocal[tp/factor,ρ[[2]]1.05 ]))}},{tp,ti factor,t factor},PlotStyle->{{Thickness[0.007],Red}}]]

plotgradv[tplotb_:0.4,endy_:ymax,color_:{Blue,,Thickness[0.005]}]:=(
tplot=tplotb/factor;
ParametricPlot[{ ((-1)/factor) Derivative[0,1][a4si][ti,(yp/factor)]/(a4si[ti,yp/factor]),v3[tplot,(yp/factor)]},{yp,0 factor,factor endy},PlotStyle->{color,{Blue},{Red,Thick},{Red}},PlotRange->{{0,1.25},{0,0.22}},ImageSize->{500,400},AspectRatio->1,PlotPoints->100,Frame->True,FrameStyle->Thickness[0.0025],FrameTicksStyle->25,FrameLabel->{Text[Style["(∂_ρ)(Subscript[ε, 0])/(Subscript[ε, 0]) ((fm^(-1)))",{FontSize->30}]],Text[Style["(Subscript[v, ρ])",{Italic,FontSize->40}]]}])

plotgradv[tplotb_:0.4,endy_:ymax,color_:{Blue,,Thickness[0.005]}]:=(
tplot=tplotb/factor;
ParametricPlot[{ ((-1)/factor) Derivative[1][a4numinti][(yp/factor)]/(a4numinti[yp/factor]-a4IRnow),v3[tplot,(yp/factor)]},{yp,0 factor,factor endy},PlotStyle->{color,{Blue},{Red,Thick},{Red}},PlotRange->{{0,1.25},{0,0.22}},ImageSize->{500,400},AspectRatio->1,PlotPoints->100,Frame->True,FrameStyle->Thickness[0.0025],FrameTicksStyle->25,FrameLabel->{Text[Style["(∂_ρ)(Subscript[ε, 0])/(Subscript[ε, 0]) ((fm^(-1)))",{FontSize->30}]],Text[Style["(Subscript[v, ρ])",{Italic,FontSize->40}]]}])

plotη:=Plot3D[#[[1]][tp/factor,yp/factor]//Abs,{tp,ti factor,t factor},{yp,0.02,factor(ymax+3)},ColorFunctionScaling->False,ColorFunction->(ColorData["Rainbow"][#3/0.4]&),PlotRange->{All,All,{0,0.4}},AxesLabel->{Text["τ(fm/c)",BaseStyle->{FontSize->30}],Text["ρ(fm)",BaseStyle->{FontSize->30,Italic}],Text[#[[2]],BaseStyle->{FontSize->30}]},LabelStyle->{Italic,FontSize->25},ImageSize->800,PlotPoints->25,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]&/@Transpose[{{((Thydroρρη[#1,#2])/(Elocal[#1,#2]))&,((Thydroθθη[#1,#2])/(Elocal[#1,#2]))&,((Thydroyyη[#1,#2])/(Elocal[#1,#2]))&},{"((Subscript[((T^visc)), ρρ])/ℰ)","((Subscript[(T^visc), θθ])/ℰ)","((Subscript[(T^visc), yy])/ℰ)"}}]

<<PhysicalConstants`;
<<Units`;
1 Giga ElectronVolt/(((Femto Meter))^3);
%/ PlanckConstantReduced;
%/ SpeedOfLight;
Convert[%,(1/(((Femto Meter))^4))]
(1/(((Femto Meter))^4))/%

plotini[color_:Blue]:=Plot[((((1/5.068)Teffa4[a4si[ti,rp/ factor]]/factor))^4),{rp,0,ymax factor},PlotRange->{{0,ymax factor},All},Frame->True,FrameStyle->Thickness[0.0025],FrameTicksStyle->25,PlotStyle->{{color,Thickness[0.005]}},ImageSize->650,FrameLabel->{Text[Style["ρ(fm)",{FontSize->40}]],Text[Style["((Subscript[T, eff])^4)((GeV^4))",{FontSize->40}]]}]

plotini[color_:Blue]:=Plot[If[rp/factor>ρ[[2]],((-27)/(5.068 8 (π^2)))a4si[ti,ρ[[2]]](factor^(-4)),((-27)/(5.068 8 (π^2)))a4si[ti,rp/ factor](factor^(-4))],{rp,0,40},PlotRange->{{0,ymax factor},All},Frame->True,FrameStyle->Thickness[0.0025],FrameTicksStyle->25,PlotStyle->{{color,Thickness[0.005]}},ImageSize->500,FrameLabel->{Text[Style["ρ (fm)",{FontSize->30}]],Text[Style["(Subscript[ε, 0])(GeV/(fm^3))",{FontSize->30}]]}]

plothydro2D:={Plot[{((Tρρ[tp,0.02]-Thydroρρ[tp,0.02])/(Elocal[tp,0.02])),((Tθθ[tp,0.02]-Thydroθθ[tp,0.02])/(Elocal[tp,0.02])),((Tyy[tp,0.02]-Thydroyy[tp,0.02])/(Elocal[tp,0.02]))},{tp,ti,t},PlotRange->{-0.15,0.15}],Plot[{((Tρρ[tp,#])/(Elocal[tp,#])),((Thydroρρ[tp,#])/(Elocal[tp,#]))},{tp,ti,t},PlotStyle->{{Blue,Thick},{Black,Dashed}}]&/@{0.02,1,3,5,10},
Plot[{Tρρ[tp,#],Thydroρρ[tp,#]},{tp,ti,t},PlotStyle->{{Blue,Thick},{Black,Dashed}}]&/@{0.02,1,3,5,10}}

plotall[color_:Blue]:={plotini[color],plotv,plota,plotf4,plotE,plotTlocal,plotT,plotPminE,plotΔP,plotη,plothydro2D}

show[data_,color_:Blue,plots_:True,useL_:20]:=(
ti=data[[3,1,1,1]];
t=data[[3,-1,1,1]];
a4IRnow=-data[[3,1,1,-1]];
changeN[35,data[[3,1]]//Length];
If[Length[data]>9,powercoord=data[[10]]];
If[Length[data]>8,
newGrid[1.01`60,data[[9]]]
,newGrid[1.01`60,useL]];
If[L<10,a4IRnow=0];

E0=300;
If[L>10,
a4iold[ρ_]:=-4E0 WSint[ρ]/3;
a4numint=a4iold/@ρ;
a4i2=Interpolation[Thread[{ρ,a4numint}][[2;;-1]],InterpolationOrder->4];
a4inf=a4numint[[1]];
a4numinti=If[#1==∞,a4inf,a4i2[#1]]&;,
a4numinti//Clear;
a4numinti[ρ_]:=-4E0 ((2/3)+(2/3)(E^(-(ρ^2)/(((20/14))^2))))/3;
];
coordtrans=Interpolation[Table[{ρp,(Sqrt[Root[-((ρp)^powercoord)+((ρp)^powercoord) #1+(L^powercoord) ((#1)^(powercoord/2))&,2]])},{ρp,0,ρ[[2]],ρ[[2]]/1500}],InterpolationOrder->4];
ymax=3L/2;
a4si=intlist[data[[3]]];
b4si=intlist[data[[4]]];
c4si=intlist[data[[5]]];
f4si=intlist[data[[6]]];
ksii=intlist[data[[7]]];
First[Timing[vlocs=Table[{#,ρ[[i]],v3[#,ρ[[i]]]},{i,Length[ρ]}]&/@Table[data[[3,i,1,1]],{i,Length[data[[3]]]}]]];
vloc=intlist[vlocs];
findfactor;
If[plots,
{findtiso,plotall[color]},
plotini[color]]
)

Do runs

Init run nucleus

sesksi="nucleus.dat";

prerun[tini_,session_,a4initial_,a4regulator_,useL_,nz2_,nρ2_]:=(
Print["Prerunning "<>session<> " with (Subscript[a, 4])(ρ)=",a4initial["ρ"],", a4IR=",a4regulator, " and (L,nz,nρ)=",{useL,nz,nρ}];
sesksi=session;
changeN[nz2,nρ2];
newGrid[1.`60,useL];
a4i[ρ_]:=a4initial[ρ];
Monitor[initKsi[tini,bi4,bi4,f4i,a4i,{600},Table[{ρ[[it]],3.2`60},{it,nρ}],False,{1.0`60,L}]//Timing,{mcon,cons,t}];
start=600-30;
stop=a4regulator 2;
changeN[nz2,nρ2];
val=Table[(E^(Log[start] (1-x)+x Log[stop])),{x,0,1,Log[1.15`60]/Log[start/stop]}]~Join~{stop};
Monitor[initKsi[tini,bi4,bi4,f4i,a4i,val,Thread[{ρ,ksiinnum}],False,{1.15`60,L}]//Timing,{mcon,cons,t}];
start=a4regulator 1.95`60;
stop=a4regulator;
val=Table[(E^(Log[start] (1-x)+x Log[stop])),{x,0,1,Log[1.05`60]/Log[start/stop]}]~Join~{stop};
Monitor[initKsi[tini,bi4,bi4,f4i,a4i,val,Thread[{ρ,ksiinnum}],True,{1.05`60,L}]//Timing,{mcon,cons,t}];
(*Thread[{ρ,ksiinnum}]>>"ksi"<>sesksi;
Thread[{ρ,ksidtinnum}]>>"ksidt"<>sesksi;
Thread[{ρ,ksidt2innum}]>>"ksidt2"<>sesksi;*)
ksidts>>"ksis"<>sesksi;
{plotConstraint,pl2/@ksidts}
);

powercoord=20;
WSint[r_]:=NIntegrate[((1/R)/(1+(E^(((Sqrt[(r^2)+(z^2)])-R)/a))))/.a->24/20/.R->31/2,{z,0,40},WorkingPrecision->60]

E0=300;
a4iold[ρ_]:=-4E0 WSint[ρ]/3
changeN[25,35];
newGrid[1.`60,18`60]

a4numint=a4iold/@ρ//Quiet;
a4i2=Interpolation[Thread[{ρ,a4numint}][[2;;-1]],InterpolationOrder->6]
a4inf=a4numint[[1]];
a4numinti=If[#1==∞,a4inf,a4i2[#1]]&

bi4[z_,ρ_]:=-(z^4)a4numinti[ρ]/8

nfilterρ2=0;
prerun[0.3`60,sesksi,a4numinti,(4/3)E0/33,18`60,35,35]

Run nucleus

sesksi

ksidts=<<("ksis"<>sesksi);
newGrid[1.01`60,L];
ksiints=Interpolation[Thread[{ρ,#}][[2;;-1]]]&/@ksidts
ksiininf=Sum[ksidts[[i,1]](((#1-ti+(10^(-60))))^(i-1))/(i-1)!,{i,Length[ksidts]}]&
ksiinρ=If[#2==∞,ksiininf[#1],Sum[ksiints[[i]][#2](((#1-ti+(10^(-60))))^(i-1))/(i-1)!,{i,Length[ksidts]}]]&

nfilterρ2=3;
nfilterz=0;
Monitor[run[0.000125`60,0.3`60,1.8`60,bi4,bi4,f4i,(a4numinti[#]-a4IRnow&),True,ksiinρ[#1,#2]&,gridsize->{1.02`60,L},execute->True,print->True,changegrid->"changeksi",onlyrun->8,computeAH->True],{cons,mcon,t}]

Init run fluctuation

sesksi="fluc.dat";

powercoord=4;

E0=300;
changeN[35,35];
newGrid[1.`60,2`60]

a4numinti//Clear;

a4numinti[ρ_]:=-4E0 ((2/3)+(2/3)(E^(-(ρ^2)/(((20/14))^2))))/3

bi4[z_,ρ_]:=-(z^4)a4numinti[ρ]/8

nfilterρ2=0;
Monitor[initKsi[3/10,bi4,bi4,f4i,a4numinti,{0,0},Thread[{ρ,Table[5/2,{nρ}]}],True,{1.02`60,L}]//Timing,{mcon,cons,t}];
ksidts>>"ksis"<>sesksi;
{plotConstraint,pl2/@ksidts}

Run fluctuation

ksidts=<<("ksis"<>sesksi);
newGrid[1.01`60,L];
ksiints=Interpolation[Thread[{ρ,#}][[2;;-1]]]&/@ksidts
ksiininf=Sum[ksidts[[i,1]](((#1-ti+(10^(-60))))^(i-1))/(i-1)!,{i,Length[ksidts]}]&
ksiinρ=If[#2==∞,ksiininf[#1],Sum[ksiints[[i]][#2](((#1-ti+(10^(-60))))^(i-1))/(i-1)!,{i,Length[ksidts]}]]&

nfilterρ2=3;
nfilterz=0;
nfilterρ=0;
a4IRnow=0;
Monitor[run[0.0001`60,0.3`60,1.7`60,bi4,bi4,f4i,(a4numinti[#]-a4IRnow&),True,ksiinρ[#1,#2]&,gridsize->{1.02`60,L},execute->True,print->True,changegrid->"changeksi",onlyrun->8,computeAH->True],{cons,mcon,t}]

Results

datanucleus=<<"resultnucleus.dat";

datafluc=<<"resultfluc.dat";

Nucleus

WSint[r_]:=NIntegrate[((1/R)/(1+(E^(((Sqrt[(r^2)+(z^2)])-R)/a))))/.a->24/20/.R->31/2,{z,0,40},WorkingPrecision->60];

plblue=show[datanucleus,Blue,False]//Quiet

plota

plothydrorho[{-.6,1}]

plotvE

e2=v3[t,#/factor](((-27)/(5.068 8 (π^2)))a4si[t,#/ factor](factor^(-4)))&;
dt=0.001;
((e2[tp+dt]-e2[tp])/dt)/.tp->0

{gradblue=plotgradv[.6,7.2/factor,{Blue,Thickness[0.006]}],
gradblue2=plotgradv[.35,7.2/factor,{Blue,Dashed,Thickness[0.006]}]}

Fluctuation

plred=show[datafluc,Red,False]

plota

plotTlocal

plotT

{gradred=plotgradv[.6,1.2/factor,{Red,Thickness[0.006]}],
gradred2=plotgradv[.35,1.2/factor,{Red,Dashed,Thickness[0.006]}]}

Combined results

Show[plblue,plred]

Show[plblue,plred,PlotRange->{{0,9},{0,1200}}]

plgrad=Show[gradblue,gradred,gradblue2,gradred2,Plot[{(0.6-ti factor) x/3,(0.35-ti factor) x/3},{x,0,2},PlotStyle->{{Black,Thickness[0.005]},{Black,Dashed,Thickness[0.005]}}],Graphics[Arrow[{{0.4,0.18},{0.15,0.18}}]]]

<<PlotLegends`

ShowLegend[plgrad,{{{Graphics[{Blue,Line[{{-1,0},{1,0}}]}],Text["Model nucleus",BaseStyle->{FontSize->16,FontFamily->"Times"}]},{Graphics[{Red,Line[{{-1,0},{1,0}}]}],Text["Model fluctuation",BaseStyle->{FontSize->16,FontFamily->"Times"}]},
{Graphics[{Black,Line[{{-1,0},{1,0}}]}],Text["Formula (10)",BaseStyle->{FontSize->16,FontFamily->"Times"}]}},LegendPosition->{-0.4,0.45},LegendSize ->{0.8,0.35}}]

ShowLegend[plgrad,{{{Graphics[{Blue,Line[{{-1,0},{1,0}}]}],Text["Model nucleus",BaseStyle->{FontSize->20,FontFamily->"Times"}]},{Graphics[{Red,Line[{{-1,0},{1,0}}]}],Text["Model fluctuation",BaseStyle->{FontSize->20,FontFamily->"Times"}]},
{Graphics[{Black,Line[{{-1,0},{1,0}}]}],Text["Formula (10)",BaseStyle->{FontSize->20,FontFamily->"Times"}]}},LegendPosition->{-0.45,0.45},LegendSize ->{1.1,0.35}}]

Movie

framev[tp_]:=Plot3D[If[(Sqrt[(xp^2)+(yp^2)+0.01])/factor>ρ[[2]],0,v3[tp,(Sqrt[(xp^2)+(yp^2)+0.01])/factor]+0.001],{xp,-ymax factor,ymax factor},{yp,-ymax factor,ymax factor},ColorFunction->Function[{x,y,z},ColorData["Rainbow"][Abs[(z/(Max[Abs[.25]]))]]],ColorFunctionScaling->False,PlotRange->{All,All,{0,0.3}},AxesLabel->{Text["x(fm)",BaseStyle->{FontSize->30}],Text["y(fm)",BaseStyle->{FontSize->30}]},PlotLabel->Text["v",BaseStyle->{Italic,FontSize->30}],LabelStyle->{Italic,FontSize->20},ImageSize->800,PlotPoints->50,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

frameE[tp_]:=Plot3D[If[((Sqrt[(xp^2)+(yp^2)])/factor)>ρ[[2]],-((27 a4si[tp,ρ[[2]]])/((5.068` 8 (π^2)) (factor^4))),-((27 a4si[tp,((Sqrt[(xp^2)+(yp^2)])/factor)])/((5.068` 8 (π^2)) (factor^4)))],{xp,-ymax factor,ymax factor},{yp,-ymax factor,ymax factor},ColorFunction->Function[{x,y,z},ColorData["Rainbow"][Abs[(z/(Max[Abs[1210]]))]]],ColorFunctionScaling->False,PlotRange->{All,All,{0,910}},AxesLabel->{Text["x(fm)",BaseStyle->{FontSize->30}],Text["y(fm)",BaseStyle->{FontSize->30}]},PlotLabel->Text[Style["(Subscript[ε, 0])(GeV/(fm^3))",{FontSize->25}]],LabelStyle->{Italic,FontSize->20},ImageSize->800,PlotPoints->30,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

$MinPrecision=0;
Monitor[First[Timing[frames=Table[{{Show[frameE[tp],Graphics3D[Text["τ = "<>ToString[SetPrecision[tp factor,2]]<>"fm/c",{3,10,800},BaseStyle->{FontSize->25}]],ImageSize->500],Show[framev[tp],ImageSize->500]}}//TableForm,{tp,ti+0.01,t,(t-ti)/50}]]],tp]

Export["nucleus.avi",frames]

framev[tp_]:=Plot3D[If[(Sqrt[(xp^2)+(yp^2)+0.01])/factor>ρ[[2]],0,v3[tp,(Sqrt[(xp^2)+(yp^2)+0.01])/factor]+0.001],{xp,-ymax factor,ymax factor},{yp,-ymax factor,ymax factor},ColorFunction->Function[{x,y,z},ColorData["Rainbow"][Abs[(z/(Max[Abs[.12]]))]]],ColorFunctionScaling->False,PlotRange->{All,All,{0,0.12}},AxesLabel->{Text["x(fm)",BaseStyle->{FontSize->30}],Text["y(fm)",BaseStyle->{FontSize->30}]},PlotLabel->Text["v",BaseStyle->{Italic,FontSize->30}],LabelStyle->{Italic,FontSize->20},ImageSize->800,PlotPoints->50,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

frameE[tp_]:=Plot3D[If[((Sqrt[(xp^2)+(yp^2)])/factor)>ρ[[2]],-((27 a4si[tp,ρ[[2]]])/((5.068` 8 (π^2)) (factor^4))),-((27 a4si[tp,((Sqrt[(xp^2)+(yp^2)])/factor)])/((5.068` 8 (π^2)) (factor^4)))],{xp,-ymax factor,ymax factor},{yp,-ymax factor,ymax factor},ColorFunction->Function[{x,y,z},ColorData["Rainbow"][Abs[(z/(Max[Abs[1210]]))]]],ColorFunctionScaling->False,PlotRange->{All,All,{0,1210}},AxesLabel->{Text["x(fm)",BaseStyle->{FontSize->30}],Text["y(fm)",BaseStyle->{FontSize->30}]},PlotLabel->Text[Style["(Subscript[ε, 0])(GeV/(fm^3))",{FontSize->25}]],LabelStyle->{Italic,FontSize->20},ImageSize->800,PlotPoints->30,FaceGridsStyle->Directive[Black,Thickness[0.001]],AxesStyle->Thickness[0.004],BoxStyle->Dashed,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},Mesh->10]

$MinPrecision=0;
Monitor[First[Timing[frames=Table[{{Show[frameE[tp],Graphics3D[Text["τ = "<>ToString[SetPrecision[tp factor,2]]<>"fm/c",{.3,1,800},BaseStyle->{FontSize->25}]],ImageSize->500],Show[framev[tp],ImageSize->500]}}//TableForm,{tp,ti+0.01,t,(t-ti)/50}]]],tp]

frames[[15]]

Export["fluctuation.avi",frames]