2021-05-17 version.
Section 1. Numerical solving of systems of LODEs via gauge transformations and implicit Runge-Kutta method

(* y->100 (original by Yi),  this notebook uses the initial cond at y->1 *)

hh[k_,n_,x_,y_]:=Integrate[t^k*Exp[-t]*HypergeometricPFQ[{},{n},t*y],{t,0,x}]; (* define the target function *)

hh1 = hh[10, 1, 1, y] (* specify k = 10, n = 1, x = 1 *)

EQ = y^4 D[f[y], {y, 4}] +(-y^4 + (2 n + 2) y^3) D[f[y], {y, 3}] + (-y^3 x + (-k-n-3) y^3 +(n^2+n) y^2) D[f[y], {y, 2}]+ ((y^3 - n y^2) x + (-n k - 2 n) y^2) D[f[y], {y, 1}] + (k + 1) y^2 x (* the differential equation satisfied by the target function *)

EQ1 = Together[(EQ /. {k -> 10, n -> 1, x->1})/y^2](* specify k = 10, n = 1, x = 1 *)

11-13 f'[y]+y f'[y]+2 f''[y]-15 y f''[y]+4 y Derivative[3][f][y]-(y^2) Derivative[3][f][y]+(y^2) Derivative[4][f][y]

(* num-ht2/hkn-misc,  Let Q be the system for EQ.  P=Q+G'/G E where E is the identity matrix and G is the Gauge trans G=Exp[-y]*y^(-(k-n+1)) *)
myrule={n->1,k->10,x->1};
 G=Exp[-y]*y^(-(k-n+1));
P={{(-y-k+n-1)/(y),1,0,0},
{0,(-y-k+n-1)/(y),1,0},
{0,0,(-y-k+n-1)/(y),1},
{((-k-1)*x)/(y^2),((-y+n)*x+n*k+2*n)/(y^2),(y*x+(k+n+3)*y-n^2-n)/(y^2),(-k-n-3)/(y)}}

(* Generate a table of values of the integral hh. It takes some time. *)
(* hh1Table=Table[{y,N[hh[k,n,x,y] /. myrule]},{y,1,40,1}] *)
hh1Table={{1,0.07810139136088579`},{2,0.1403534462777064`},{3,0.22982285915267087`},{4,0.3546884049781732`},{5,0.5250013288862653`},{6,0.7530317625677715`},{7,1.0536689374445893`},{8,1.4448823960377706`},{9,1.9482522503743314`},{10,2.58957747265436`},{11,3.399572237193556`},{12,4.414661473030726`},{13,5.67788804349795`},{14,7.2399453532959654`},{15,9.160350706917452`},{16,11.508776417322787`},{17,14.366557504368297`},{18,17.82839684353342`},{19,22.00429084314467`},{20,27.021701160033906`},{21,33.02800062830995`},{22,40.19322449412041`},{23,48.713161243043935`},{24,58.812820799972776`},{25,70.75032169981955`},{26,84.82124299897478`},{27,101.36349125221419`},{28,120.76273785013576`},{29,143.45848743319286`},{30,169.95084400766547`},{31,200.80804782712895`},{32,236.67486311391784`},{33,278.2819043258564`},{34,326.455996974932`},{35,382.13167803117784`},{36,446.36395075563235`},{37,520.342419463968`},{38,605.4069412952418`},{39,703.064944621151`},{40,815.0105773587113`}};

hh1GTable=Table[{hh1Table[[i]][[1]],hh1Table[[i]][[2]]*N[G /. myrule /. {y->hh1Table[[i]][[1]]}]},{i,1,Length[hh1Table]}];
hh1func=Interpolation[hh1Table]; hh1Gfunc = Interpolation[hh1GTable];

(* For[i=1, i≤Length[hh1Table],i++,Print[hh1Table[[i]][[1]],"  ",hh1Table[[i]][[2]]]] *)

init1 = Table[N[D[hh1, {y, i}] /. y -> 1, $MachinePrecision], {i, 0, 3}] (* compute the initial value vector at y = 1 *)

init1b=N[G /. myrule /. {y->1}]*init1

g[y]= {g1[y],g2[y],g3[y],g4[y]};

dg=D[g[y],y];
EQ1b=P.g[y];

sys=Join[Table[dg[[i]]==EQ1b[[i]],{i,1,4}] /. myrule,{g1[y] == init1b[[1]],g2[y]==init1b[[2]],g3[y]==init1b[[3]],g4[y]==init1b[[4]]}/. {y->1}]

solb=NDSolve[sys ,{g1,g2,g3,g4},{y,1,150},Method->{"FixedStep", Method->{"ImplicitRungeKutta","DifferenceOrder"->10}}, StartingStepSize->1/10]

(g1[100] /. solb) /N[G /. myrule /. {y->100},100]

N[hh1 /. y-> 100]

N[hh1 /. y-> 100] (* this is the correct value of hh1 at y = 100, which is different from that of numerical solving of the ode *)

(g1[119] /. solb) /N[G /. myrule /. {y->119},100]

N[hh1 /. y-> 119] (* by gauge transformation, solb has one digit 
precision as hh1 in [100, 119] *)

solb=NDSolve[sys ,{g1,g2,g3,g4},{y,1,150},Method->{"FixedStep", Method->{"ImplicitRungeKutta","DifferenceOrder"->10}}, StartingStepSize->1/10^3] (*decreasing StartingStepSize to 10^(-3), the computational time becomes longer. However,  decreasing StartingStepSize to 10^(-4) does not help. Change FixedStep to DoubleStep does not help. Change DifferenceOrder also does not help. *)

(g1[140] /. solb) /N[G /. myrule /. {y->140},100]

N[hh1 /. y-> 140](* by gauge transformation and decreasing StartingStepSize to 10^(-3), solb has one digit 
precision as hh1 in [100, 140] *)

2. Explicit RungeKutta

solc=NDSolve[sys ,{g1,g2,g3,g4},{y,1,40},Method->{"TimeIntegration"->"ExplicitRungeKutta"},StartingStepSize->1/10^3]

Plot[{Log[(g1[y] /. solb)/N[G /. myrule ]],Log[hh1func[y]],Log[(g1[y] /. solc)/N[G /. myrule]]},{y,1,35}]

(* Plot[((g1[y] /. solb)/N[G /. myrule])/hh1-1,{y,1,35}] *)

3. IDA (DFE?) multistep method

sold=NDSolve[sys ,{g1,g2,g3,g4},{y,1,40},Method->{IDA}]

Plot[{Log[(g1[y] /. sold)/N[G /. myrule ]],Log[hh1func[y]]},{y,1,35}]

4. Solving the boundary value problem

(* Solving the boundary value problem for Gauge transformed system. *)
Clear[init];
init[y_]:= N[hh[k,n,x,y] /. myrule, $MachinePrecision] *N[G /. myrule, $MachinePrecision]
boundaryCondition=Table[g1[y] == init[y],{y,{0.98,0.99,10.1,10.01}}]
sysBoundary=Join[Table[dg[[i]]==EQ1b[[i]],{i,1,4}] /. myrule,boundaryCondition]

sole=NDSolve[sysBoundary ,{g1,g2,g3,g4},{y,1,9.99}]   (* Change in 2021-05-17 *)

boundaryCondition

Plot[{Log[(g1[y] /. sole)/N[G /. myrule ]],Log[hh1func[y]]},{y,1,10}]

(* solving the boundary problem for ODE with no Gauge transformation, Yaguchi case 2021-05-17 *)
finit[y_]:= N[hh[k,n,x,y] /. myrule, $MachinePrecision]
fboundaryConditionY=Table[f[y] == finit[y],{y,{0.98,0.99,10.0,10.01}}]
fsysBoundaryY=Join[{EQ1==0}/. myrule,fboundaryConditionY]
(* save in num-ht3/Figs/2021-05-17-math-boundary.eps *)

(* solfY=NDSolve[fsysBoundaryY ,f,{y,1,9.99},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement"}}}] *)
 (* solfY=NDSolve[fsysBoundaryY ,f,{y,1,9.99},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MinPoints"->1000}}}] *)
solfY=NDSolve[fsysBoundaryY ,f,{y,1,9.99}]
Plot[{f[y] /. solfY,hh1func[y]},{y,1,9.99}]

(* solving the boundary problem for ODE with no Gauge transformation *)
finit[y_]:= N[hh[k,n,x,y] /. myrule, $MachinePrecision]
fboundaryCondition=Table[f[y] == finit[y],{y,{1,3,6,10}}]
fsysBoundary=Join[{EQ1==0}/. myrule,fboundaryCondition]
solf=NDSolve[fsysBoundary ,f,{y,1,9.99}]  (* 2021-05-17 version *)

Plot[{Log[(f[y] /. solf)],Log[hh1func[y]]},{y,1,10}]

(* solving the boundary problem for ODE with no Gauge transformation 3 *)

fboundaryCondition3=Table[f[y] == finit[y],{y,{4,6,8,10}}]
fsysBoundary3=Join[{EQ1==0}/. myrule,fboundaryCondition3]
solf3=NDSolve[fsysBoundary3 ,f,{y,4,10}]
Plot[{Log[(f[y] /. solf3)],Log[hh1func[y]]},{y,4,10}]

5. misc

hs=h0+h1*y+h2*y^2/2 + h3*y^3/6+h4*y^4/24

Expand[(hs /. {y->2*y}) - 2*hs + 2*(hs /. {y->-y}) - (hs /. {y->-2*y})]

Expand[(hs /. {y->2*y})-4*hs+6*h0-4*(hs /. {y->-y})+(hs /. {y->-2*y})]

tmp1=Table[N[finit[y]],{y,{3.98,3.99,30.0,30.01,40,40.01}}];
tmp1

ToString[NumberForm[tmp1,15]]