eval(lhs(eq),{sin=0,cos=1});
eval(lhs(eq),{sin=1,cos=0});

Consider these 4 procedures, the first 3 with seq, the 4th with $ :

p1:=proc(n) local j,a,b; a:=j; b:=seq(a,j=1..n); b end proc;
p1(3); # j,j,j

p2:=proc(n) local j,a,b; a:=j; b:='seq'(a,j=1..n); b end proc;
p2(3); #Unevaluated seq.
%; #now fully evaluated

p3:=proc(n) local j,a,b; a:=j; b:=seq(eval(a),j=1..n); b end proc;
p3(3); # 1,2,3
## And now using $ :
q:=proc(n) local j,a,b; a:=j; b:=a$j=1..n; b end proc;
q(3); # 1,2,3

## You can replace seq with add or mul to get the same behavior in p1,p2,p3.
## You can replace $ with sum or product in q to get the same behavior.

I think this seq behavior is a combination of evaluation rules in procedures and the special evaluation rules for seq (similarly for add and mul).
### Try in p3 the following variation, where evaluation is only done to 1 level:
p3:=proc(n) local j,a,b; a:=j; b:=seq(eval(a,1),j=1..n); b end proc;
p3(3); #now j,j,j
Try replacing eval(a,1) by eval(a,2) and you have 1,2,3.
##########################
Now try at the interactive level:
restart;
a:=j;
eval(a,1); #Just to see
seq(eval(a,1),j=1..3); #j,j,j
## This latter behavior mimics the behavior of seq in the procedure p1.
###
### In the help page for seq you can find a description of how seq works, where it is constructed as a loop.
Adapting that to the situation in a procedure, where evaluation is done to one level only, we get for our present example:
restart;
a:=j;
S:=NULL;
old:=j;
for j from 1 to 3 do S:=eval(S,1),eval(a,1) end do;
j:=eval(old,1);
S;  #j,j,j

It may help to compare to sin^2 as done here:

f:=x^2;
f(3);
g:=sin^2;
g(3);
### It may also be illuminating to try this:
x^2(3);

The explanation for the output x^2 is that 2(3) is computed first. Numbers can act as constant functions, thus 2(3) = 2.

To get around the latter example use parentheses:
(x^2)(3);

We first show that the expression u below is independent of the variables:
u:=-Omega*a*sqrt(2)*sqrt(-Omega^2*a^2-2*k*m+sqrt(Omega^2*a^2*(Omega^2*a^2+4*k*m)))/(-Omega^2*a^2+sqrt(Omega^2*a^2*(Omega^2*a^2+4*k*m)));

## Clearly u only depends on the products Omega*a and k*m, so we need only consider:
w:=subs(Omega=1,m=1,u);
simplify(diff(w,a));
simplify(diff(w,k));
# Thus w is independent of a and k.
# Thus its value is:
simplify(eval(w,{a=1,k=1}));

foo := proc({ctx :: list := []}) bar(ctx = ctx) end proc:
bar := proc({ctx :: list := []}) nops(ctx) end proc:

When foo is called as in
foo(ctx = [a,b,c]);

then ctx is actually inside the procedure foo given the value [a,b,c] all over.
That means that inside foo you now have bar([a,b,c] = [a,b,c]).
Thus bar is not called with a list as an argument, so the value of nops([]) is returned, i.e. 0.

This is most definitely the intended design!

See
?Procedure Parameter Declarations

where in the section Keyword Parameters you find:

A keyword parameter receives a value when an argument of the form keyword=value appears in a function call. The left hand side of the argument specifies the keyword parameter that is to receive a value, and the right hand side specifies the value it is to receive.

I made several changes to your worksheet:
(1) Digits was misspelled. I set it at 15.
(2) I used the option parameter, so dependence of omega could be examined in an animation. After all you have omega=10^6 (!!)
(3) I use maxfun=0 (which means maxfun=infinity).
(4) No range argument.
(5) I plot with odeplot.
##
restart;
Digits:=15:
x(0):=-1:y(0):=1:z(0):=conjugate(y(0)):N:=10:Delta:=5: N1:=1+2*N:M:=sqrt(N*(N+1)):
ini1:=u(0)=Re(y(0)), v(0)=Im(z(0)),w(0)=x(0);
#omega:=10^(6): #NOT assigned
dsys1 :=diff(w(t),t)=-(N1+M*cos(2*omega*t))*w(t)-1+2*u(t)*cos(2*omega*t)+2*v(t)*sin(2*omega*t), diff(u(t),t)=-N1*u(t)+Delta*v(t)-2*M+(2*M*u(t)-N1-w(t))*cos(2*omega*t)-2*M*v(t)*sin(2*omega*t), diff(v(t),t)=-N1*v(t)-Delta*u(t)-2*M+(2*M*u(t)-N1-w(t))*sin(2*omega*t)+2*M*v(t)*cos(2*omega*t);

dsol1 :=dsolve({dsys1,ini1},numeric, parameters=[omega],abserr=1e-9, relerr=1e-8,maxfun=0);
## The procedure p is used for animation:
p:=proc(omega) dsol1(parameters=[omega]);
   plots:-odeplot(dsol1,[[t,u(t)],[t,v(t)],[t,w(t)]],0..1,numpoints=10000)
end proc;
plots:-animate(p,[omega],omega=1000..5000,size=[1800,default]); #Adjust size to your screen size
##Notice that the wiggles are almost gone at omega=5000.
dsol1(parameters=[10^6]); #Setting omega to 10^6.
plots:-odeplot(dsol1,[[t,u(t)],[t,v(t)],[t,w(t)]],0..1,size=[1800,default],numpoints=10000);
##No perceptible wiggles.

You could do:

{(x-2)/3 , (y-1)/4 , (z-3)/3}=~t;
solve(%,{x,y,z});
L:=subs(%,[x,y,z]);
plots:-spacecurve(L,t=-5..5,thickness=3,color=red);

Add unevaluation quotes:
g := unapply('unapply(theta*x,x)',theta);

Without the quotes the inner unapply is done first resulting in a procedure x->theta*x. The next unapply won't see the theta.
It would be like doing g:=unapply(x->theta*x,theta);
In doing this unapply will evaluate (x->theta*x), but nothing is done, quite as it should be!

What to do after a singularity (a genuine one) happens is up to you it seems to me.
You could try this, where the event that T(t)-0.001=0 triggers the change T(t)=.001 and n(t)=.001.
That is clearly arbitrary. You can experiment with other values for the changes.

resE:=dsolve(sys_ODE1, numeric,events=[[T(t)-.001,[T(t)=.001,n(t)=.001]]]);
plots:-odeplot(resE,[[t, log10(n(t))], [t, log10(T(t))]],0..50,size=[1000,default]);

Try
restart;
example1 := Matrix([[1, 2, 3], [4, 5, 6], [6, 7, 8]]);
testproc := proc (A) whattype(A) end proc; #No print
testproc(example1);
whattype(example1);
print(%);
## The explanation is that Matrix itself is a procedure (used for constructing matrices):
type(Matrix,procedure);
# If you want to see it:
showstat(Matrix);

It should be
plots[animate](plot, [[f(x), h(x), x = 0 .. 1]], h0 = 0 .. 1);
as in the example:
plots[animate](plot, [[x^2+h0, sin(h0*x), x = 0 .. 1]], h0 = 0 .. 1);

Here is a solution, where I set infinity to 5.
restart;
ode:=diff(f(x),x,x,x)+f(x)*diff(f(x),x,x)+1-diff(f(x),x)^2=0;
bcs:=f(0)=0,D(f)(0)=0,D(f)(infinity)=1;
res:=dsolve(eval({ode,bcs},infinity=5),numeric);
plots:-odeplot(res,[seq([x,diff(f(x),[x$i])],i=0..2)],0..5);

#I'll leave the saving to somebody else.
############################# Elaborating and using shooting too ###################
restart;
Digits:=12: #fsolve problems at Digits=15
ode:=diff(f(x),x,x,x)+f(x)*diff(f(x),x,x)+1-diff(f(x),x)^2=0;
bcs:=f(0)=0,D(f)(0)=0,D(f)(infinity)=1;
res:=dsolve(eval({ode,bcs},infinity=5),numeric);
plots:-odeplot(res,[seq([x,diff(f(x),[x$i])],i=0..2)],0..5,color=[red,blue,green]); p5:=%:
res(0);
f2_res:=subs(res(0),diff(f(x),x,x));
##Nothing new so far.
#############################
##Now shooting with f''(0)=f2 as shooting parameter.
## We will determine f2 s.t. f'(inf)=1.
ics:=f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=f2; #Initial conditions.
respar:=dsolve({ode,ics},numeric,parameters=[f2],output=listprocedure);
F1:=subs(respar,diff(f(x),x)); #Procedure to compute f'(x).
inf:=10: #Infinity set at 10.
## Procedure Q enabling fsolve to find f2:
Q:=proc(f2) respar(parameters=[f2]); F1(inf)-1 end proc;
plot(Q,1.1..1.3,-1..1); #We see two zeros, i.e. solutions for f2!
plot(Q,1.14..1.16,-1..1); #Zooming in.
f2_1:=fsolve(Q,1.15..1.152);
respar(parameters=[f2_1]);
plots:-odeplot(respar,[seq([x,diff(f(x),[x$i])],i=0..2)],0..inf); p1:=%: #Not the one we want.
f2_2:=fsolve(Q,1.23);
respar(parameters=[f2_2]);
plots:-odeplot(respar,[seq([x,diff(f(x),[x$i])],i=0..2)],0..inf); p2:=%: #Looks good!
plots:-display(p1,p2); #Two different solutions on x=0..inf.
plots:-display(p2,p5,view=[0..5,0..5]); #Agreement between p2 and the bvp solution p5.
## Looking at the dependence of the two zeros on inf:
Qanim:=proc(f2,INF)
   if not type([_passed],list(numeric)) then return 'procname(_passed)' end if;
   respar(parameters=[f2]);
   F1(INF)-1
end proc;
plots:-animate(plot,[Qanim(f2,INF),f2=0.8..1.3,-2..2],INF=5..30,frames=10,trace=10);

The remarkable thing here is the fact that the second zero seems independent of the value of inf!

op combined with convert D is not all that bad:
expr := diff(f(x,y),x,y$2);
E:=convert(expr,D);
op(E); #variables
op([0,1],E); #function
op(op([0,0],E)); #Derivatives
###
Since you are already aware of PDEtools:-difforder I cannot suggest anything else.
indets(expr,name); #variables
indets(expr,function(name)); #function
PDEtools:-difforder(expr); #total order
PDEtools:-difforder(expr,x); #order wrt. x
PDEtools:-difforder(expr,y); #order wrt. y
###
Procedure for the first:
pp:=proc(u) local E,var,f,L;
    E:=convert(u,D);
    var:=op(E);
    f:=op([0,1],E);
    L:=op(op([0,0],E));
  f(var),[seq(var[i],i=L)]
end proc;

a:= b(t)+diff(b(t),t)+diff(b(t),t$2)+diff(b(t),t$3);
convert(a,D);
subs(D(b)(t)=b_symbol_diff,%);
convert(%,diff);
#########################
Another idea:

S,R:=selectremove(type,a,specfunc(Not(specfunc(diff)),diff));
R+subs(diff(b(t),t)=b_symbol_diff,S);

You are being told quite explicitly what the problem is: Division by zero.
Indeed you have this line
bcs2 := eval[recurse](convert(ds[1], D), `union`({y = 0}, bcs3));
If you look at
convert(ds[1], D);
you will see divisions by various positive powers of y. You intend to evaluate that at y=0. Clearly a problem you cannot run away from.

(And what is the value of L?)