Using NeumannValue in finite element method to replace a derivative boundary condition in NDSolve to overcome a computational overflow

Question

I am trying to implement and evaluate the function GeneratePsiOfX. The definition of GeneratePsiOfX is:

GeneratePsiOfX[psi0_?NumericQ,xmin_?NumericQ,xmax_?NumericQ,maxsteps_?NumericQ,method_]:=Module[{ode,ID,ORD,eq1,ss,as,Eqs,Soln,Soln2,tmp},
ode=2 (\[Psi]')[x]+x (\[Psi]'')[x]+3x E^\[Psi][x] Sqrt[\[Pi]] Erf[Sqrt[\[Psi][x]]]-2 Sqrt[\[Psi][x]] (3+2 \[Psi][x])x;
ID={\[Psi]0->psi0,\[Psi]p0->0};
ORD=10;
eq1=\[Psi]->Function[x,Sum[a[i]x^i,{i,0,ORD}]/.{a[0]->\[Psi]0,a[1]->\[Psi]p0}/.ID];
ss=Series[ode/.eq1,{x,0,ORD}];
(*Print[ss];*)
as=Table[a[i],{i,2,ORD}]; 
Eqs=Table[SeriesCoefficient[ss,i]==0,{i,1,ORD}];
(*Print[Eqs];*)
seriesSoln=Solve[Eqs,as]//Flatten;
Soln=\[Psi][x]/.eq1/.seriesSoln//N;
ode1=2 (\[Psi]')[x]+x (\[Psi]'')[x]+3x E^\[Psi][x] Sqrt[\[Pi]] Erf[Sqrt[\[Psi][x]]]-2 Sqrt[\[Psi][x]] (3+2 \[Psi][x])x==0;
xstart=xmin;
xstop=xmax;
Soln2=NDSolve[{ode1,\[Psi][xstart]==Soln/.x->xstart,\[Psi]'[xstart]==D[Soln,x]/.x->xstart},\[Psi][x],{x,xstart,xstop},MaxSteps->maxsteps,Method->method][[1,1,2]];
Piecewise[{{Soln,x<=xstart},{Soln2,x>xstart}}]
]

Upon evaluating GeneratePsiOfX as follows:

GeneratePsiOfX[15,0.001,4,10^7,{"StiffnessSwitching"}]

I get the following error:

General::ovfl: Overflow occurred in computation.

So, I decided to switch to Finite Element Method instead of Stiffness Switching as below:

GeneratePsiOfX[15,0.001,4,10^7,{"FiniteElement"}]

I now get the following error:

NDSolve::fembderiv: The expression (\[Psi]')[0.001]==-29568.3 given as a spatial boundary condition for the possibly automatically chosen finite element method should not have explicit derivatives of the dependent variables. NeumannValue should be used to specify spatial derivatives at the boundary.

I read the official documentation of NeumannValue but it is a little difficult to understand. I tried reading a few posts on this website as well. Here is what I think is required. I modified the NDSolve function within the GeneratePsiOfX function to include NeumannValue as follows:

Soln2 = NDSolve[{ode1, \[Psi][xstart] == Soln /. x -> xstart, NeumannValue[D[Soln,x]/.x->xstart,x==xstart]}, \[Psi][x], {x, xstart, xstop}, MaxSteps -> maxsteps, 
    Method -> method][[1, 1, 2]];

However, when I implement the following:

GeneratePsiOfX[15, 0.001, 4, 10^7, {"FiniteElement"}]

I get the error:

NDSolve::deqn: Equation or list of equations expected instead of NeumannValue[-29568.3,x==0.001] in the first argument {3 E^\[Psi][x] Sqrt[\[Pi]] x Erf[Sqrt[\[Psi][x]]]-2 x Sqrt[\[Psi][x]] (3+2 \[Psi][x])+2 (\[Psi]')[x]+x (\[Psi]'')[x]==0,\[Psi][0.001]==10.8068,NeumannValue[-29568.3,x==0.001]}.

So, Mathematica first tells me to use NeumannValue but then when I do so, it tells me that's not what it expects. I am a little confused as what to do. Also, I am using Mathematica v.13.0.1.0.