You can use the parameters option.
Leave out the assignment to n.
Then replace the dsolve statement with:
sol := dsolve({gl1, gl2, gl3, gl4, init1, init2, init3, init4}, numeric,parameters=[n]);
### Example of setting n:
sol(parameters=[4000]);
sol(.5);
odeplot(sol, [x(t), y(t)], t = 0 .. 6.7);
### Now for your actual question. Use this simple procedure:
psol:=proc(n) sol(parameters=[n]); sol end proc;
## Then display the curves for n=1500..9000 (in steps of 1500):
display(seq(odeplot(psol(nn),[x(t), y(t)], t = 0 .. 6.7),nn=1500..9000,1500));

A quote from the help page for Numerical Summation:
   "In the case of infinite sums, Levin's u-transform is used."

You can see the procedure used:

showstat(`evalf/Sum/LevinU_Symb/1_inf`);
showstat(`evalf/Sum/LevinU`);

You have a space after PolynomialInterpolation and before ([Points], y). That space is interpreted by the 2D parser as multiplication, which obviously makes no sense.

The matrix options are unchanged by elementwise operations.
Try:
MatrixOptions(DC);
      shape = [], datatype = complex[8], storage = rectangular, order = Fortran_order
MatrixOptions(DX);
      shape = [], datatype = complex[8], storage = rectangular, order = Fortran_order
# Notice in particular that the datatype is complex[8] in both cases.
## Same result for
MatrixOptions(Re~(DX));
## You can give your matrix explicitly the datatype float:
DX1:=Matrix(abs~(DC),datatype=float);
MatrixOptions(DX1);
       shape = [], datatype = float[8], storage = rectangular, order = Fortran_order
## Finally, by using simplify/zero as in Adri van der Meer's answer, you get:
MatrixOptions(simplify(DX,zero)); #or DC, same answer for the options
       shape = [], datatype = anything, storage = rectangular, order = Fortran_order

I shall only deal with the first part, that of applying a matrix of operators Q to N.
You could use the second version of Q containing X,Y,Z instead of the differential operators. Then find B:=Q.N;

Now use these 3 procedures px,py,pz:
px:=proc(u) global X,x;
  if not type(u,`*`) then return u end if;
  if member(X,{op(u)}) then diff(u/X,x) else u end if
end proc;
py:=subs(X=Y,x=y,eval(px));
pz:=subs(X=Z,x=z,eval(px));
##
Then use map:
R:=map(px@py@pz,B);
## Now you have your matrix R.

1. You have D(phi) = 1 in your boundary condition. You probably meant either D(phi)(0) = 1 or D(phi)(1) = 1.
2. What is the value of n and upsilon in approxsoln?
3. What are the values of h and lambda (assuming that omega2 is to be determined)?

Your problem is of the form

eq:=a*y^beta+b*y+c=0;

where beta is just known to be nonnegative. You would certainly have luck with beta = 0, 1, 2, 3, and 4 and for those values of beta for which the equation can be reformulated as a polynomial equation of degree up to 4 in some variable z. As an example of the latter take beta = -1/2 and use z=y^(1/2).
For general beta >=0 I don't see any way of getting a formula for a solution.

If there is a singularity, you obviously cannot get rid of it. The error message says "probably a singularity" though.
But it is not unlikely that there actually is one.
Try first with less problematic numbers:
test1:=evalindets(test,float,1);
sol1 := dsolve({test1, H(0) = 1}, numeric);
plots:-odeplot(sol1,[z,H(z)],0..0.32);
##The right hand side of test1 is larger than the right hand side of this:
test2:=diff(H(z), z) = H(z)^2;
sol2 := dsolve({test2, H(0) = 1});
## which has the solution H(z) = -1/(z-1), i.e. a singularity at z=1. Thus test1 with H(0)=1 must have a singularity before z=1.
######
## Now look at test and likewise compare it to test0:
test0:=diff(H(z), z) = 6.534101519*10^17*H(z)^2;
sol0:=dsolve({test0, H(0) = 2.268308490*10^(-18)});
solve(denom(rhs(sol0))=0,z);
evalf(%); # Answer 0.6747020100
# We conclude that test with H(0)=2.268308490*10^(-18) must have a singularity before z = 0.6747020100
## Please notice that test2 and test0 both were solved using the symbolic solver and that those equations easily could have been solved by hand. Thus there is no doubt remaining regarding the existence of the singularity: It is not a numerical problem!
##If you want to push the upper bound for the singularity further down you could compare it to:
test01:=diff(H(z), z) = 6.534101519*10^17*H(z)^2*(1.+z)^(5/2);
sol01:=dsolve({test01, H(0) = 2.268308490*10^(-18)});
## Here you will see that the singularity will be below z = .413957797.

There is ListTools:-PartialSums which might be useful. It provides all partial sums of a list.

Example:
L:=[seq(a[k],k=1..10)];
ListTools:-PartialSums(L);
#Procedure:
PS:=proc(an::algebraic,r::name=range(integer)) local n,N1,N2,k;
  n:=lhs(r);
  N1,N2:=op(rhs(r));
  ListTools:-PartialSums([seq(eval(an,n=k),k=N1..N2)])
end proc;
 
PS(n^2,n=2..6);
PS(a[i],i=0..3);

Use combine on all three as in:

exp(varepsilon[10])*a[2]*varepsilon[6]*varepsilon[9]^2/exp(varepsilon[3]);
combine(%);

Begin by setting Digits:=30;
Then proceed as in your worksheet with defining all quantities up to and including eqX.

Then do
eqX:=fnormal(eqX,25);
x:=Vector([x1, x2, x3, x4, x5, x6](t));
dxdiff:=diff~(x,t);
##You cannot impose conditions on the derivatives of the xi's. So use:
ics := {x1(0) = 10^(-6), x2(0) = 10^(-6), x3(0) = 10^(-6), x4(0) = 10^(-6), x5(0) = 10^(-6), x6(0) = 10^(-6)};
sys:=convert(Equate(dxdiff,eqX.x),set);
## First a numerical solution:
res:=dsolve(sys union ics,numeric,abserr=1e-15,relerr=1e-15);
plots:-odeplot(res,[t,x1(t)],0..0.1);
##Next a formula for the solution:
sol:=evalf(dsolve(sys union ics,convert(x,list),method=laplace));
plot(subs(sol,x1(t)),t=0..0.1);

MaplePrimes16-04-21ricatti.mw

Using only directly available Maple code you could do:

u:=polar(5,30/180*Pi)+polar(4,60/180*Pi);
evalc(u);
convert(%,polar);
simplify(%);
evalf(%);
##Of course the answer uses radians in the second argument.

You can use evalindets:
a:= [2+I,2,5];
evalindets(a,And(complex,Not(realcons)),0);
## And for your example in the file:
evalindets(SolPointtheta[1],And(complex,Not(realcons)),0);

For your particular example you can get an exact answer and use spacecurve on that. But you can also use dsolve/numeric and odeplot:

restart;
EF := {2*(diff(w[2](t), t)) = 10, diff(w[1](t), t) = sqrt(2/w[1](t)), diff(w[3](t), t) = 0};
sol:=dsolve(EF union {w[1](0) = 1, w[2](0) = 0, w[3](0) = 0}) ; #Exact
DEtools[remove_RootOf](sol[1]); #To see the RootOf equation
allvalues(sol); #3 solutions, 2 imaginary
sol1:=op(remove(hastype,[%],imaginary));
sol2:=subs(sol1,[w[1](t),w[2](t),w[3](t)]);
plots:-spacecurve(sol2,t=0..10,coords=spherical); #Using spacecurve
L:=[r*cos(theta)*sin(phi),r*sin(theta)*sin(phi),r*cos(phi)]; #x,y,z in spherical coords
plots:-spacecurve(subs({r=sol2[1],theta=sol2[2],phi=sol2[3]},L),t=0..10); #Just to confirm the result above
## Now numerical solution:
res:=dsolve(EF union {w[1](0) = 1, w[2](0) = 0, w[3](0) = 0},numeric) ;
plots:-odeplot(res,subs({r=w[1](t),theta=w[2](t),phi=w[3](t)},L),0..10);

You could use mtaylor and freeze as follows.
Let expr be only the left hand side of your big expression.
Borrowing some symbols from acer do:

L:={a,b,c,u,v,w,psi,phi,theta,varsigma,tau,upsilon}(t):
LL:=L union map(diff,L,t) union map(diff,L,t,t):
S:=LL =~freeze~(LL);
mtaylor(subs(S,expr),rhs~(S),2);
thaw(%);