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?)

restart;
sys:= diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))+S(t)*(-z*eta*alpha*K(t)^2+(-z*eta*alpha*S(t)-(eta*alpha^2*S(t)^2-2*N*C[max]*w*eta*alpha*K(t)+((-N*w+z)*alpha+N*C[max]^2*w*eta)*w*N)*upsilon)*K(t)+N*S(t)*w*alpha*upsilon*(N*w-z))/((K(t)^2*alpha*z+3*K(t)*S(t)*alpha*z+(2*S(t)*z*alpha+upsilon)*S(t))*w*N);
eq1:= -c2 + (K(t2)*S(t2)+w*N*S(t2))*z=0;
ics:=S(0)=100,K(0)=20;
params:=indets({sys,eq1},name) minus {t,t2};
##Now I set all parameters to 1 since I don't know what they are:
vals:=params=~1;
##Use your own of course.
eval([sys,eq1],vals);
eq2:=subs(t2=t,%[-1]);
sys2:=%%[1..2][];
#Because of my parameters eq2 never becomes satisfied so for illustration I used the equation
# lhs(eq2)=200 instead
res:=dsolve({sys2,ics,T(0)=0},numeric,discrete_variables=[T(t)::float],events=[[lhs(eq2)=200,T(t)=t]]);
res(5);
plots:-odeplot(res,[t,lhs(eq2)-200],0..0.1);

In Maple go to the menu Tools:  Tools/Tutors/Linear Algebra/Gaussian Elimination ...

Here is another version. It uses 'satisfies'.

createModule2 := proc(A::Matrix(square))
    local dim;
    dim := LinearAlgebra:-RowDimension(A);
    module()
        export det;
        det := (x::And(Matrix(square),satisfies(x->op([1,1],x)=dim))) -> LinearAlgebra:-Determinant(x)
    end module
end proc:

I shall just call your attention to two procedures CurveFitting:-PolynomialInterpolation and Statistics:-Fit.

data:=[[1030, 0], [380, 34], [270, 73], [240, 150], [85, 700], [22, 2000], [12, 5000]];
##This one uses collocation, i.e, the curve passes through all the points.
p:=CurveFitting:-PolynomialInterpolation(data, x);
plots:-pointplot(data,symbolsize=20); p1:=%:
plots:-display(plot(p,x=0..1030),p1); # Not great
## Best fit to exponential:
q:=Statistics:-Fit(a*exp(-b*x),data,x);
plots:-display(plot(q,x=0..1030),p1); #Not particularly great either

So think of some reasonable model and use Statistics:-Fit.