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(%);

I don't believe that there is any way of getting a formula for the general solution, nor is there a way to find a formula for the ivp {eq, y(0) = 283.8}.

But you can easily solve numerically:

res:= dsolve({eq, y(0) = 283.8},numeric);
plots:-odeplot(res,[x,y(x)],-1500..1465);


####
# Note:
dsolve(eq);
returned nothing (NULL) in Maple 17.02 as well as in Maple 2016 on my computer.

Use fnormal.
See the help page ?fnormal.
But here is a simple example (Digits is default 10 on my computer):
restart;
Digits; #10
u:=a*1e-3+b*1e-9+c*1e-15;
fnormal(u);
fnormal(u,11);
fnormal(u,17);

The D you get in the following is the usual D.
restart;
b:=<b1(t),b2(t)>;
f(b); #Not a container (relevant for elementwise ~)
##The only container in the following is b:
dfdb:=Physics[diff]~(f(b),b); #Output a vector with two identical elements! Is that intended?
dfdb[1];
               D(f)(<b1(t), b2(t)>)  ##Output (same as dfdb[2])

If we now define f as
f:=b->b[1]^2+b[2];

then notice that even now the following doesn't evaluate to anything other than the input:
D(f)(b);
               D(f)(<b1(t), b2(t)>)  ##Output
###
If you make f a function of two variables instead of a function of a vector you get more of what you want:
restart;
b:=<b1(t),b2(t)>;
dfdb:=Physics[diff]~(f(b1(t),b2(t)),b);
f:=(b1,b2)->b1^2+b2;
dfdb;

I think the explanation is that the harm has already been done before evalf gets to the number:
restart;
5^1.25; #10 digits with default settings of Digits
'5^1.25'; #Unevaluation quotes don't help: automatic simplification.
5^(5/4); #Still symbolic

So raise Digits:
Digits:=30:
evalf(5^1.25);