Even for mm=1 and nn=1 it is impossible:

f:=t->add(A[m]*cos(m*w*t+m*t0),m=0..1);
g:=phi->add(B[n]*cos(n*phi+n*phi0),n=0..1);

fg:=unapply(combine(f(t)*g(phi)),t,phi);
h:=(t,phi)->add(add(C[m,n]*cos(m*w*t+n*phi+cmn),n=0..1),m=0..1);
eq:=h(t,phi)-fg(t,phi)=0;
diff(eq,phi);
diff(%,t);
subs(w*t=s-phi,%);
diff(%,phi);
#Result: A[1]*B[1]*sin(-s+2*phi-t0+phi0)*w = 0, which means that one of A[1], B[1], w must be zero.

You may want to take a look at ImportMatrix as it has an option 'skiplines=n'.

So if your file looks like this:

This is text
bla bla bla
ååååøøøææøælopk567
----------------
78 89 45
56 7 -98

and is saved in F:/test.txt  then do

ImportMatrix("F:/test.txt",delimiter=" ",skiplines=4);

There are a couple of obvious problems:

(1) Presumably you want to use the package T, so you need to do with(T); or use S:=T:-sin(t);

(2) You wrote DS:=diff(DS,t): This would make DS=0 and subsequently DDS=0. You probably meant DS:=diff(S,t):

However, the inhomogeneous term is irrelevant when it comes to the angular eigenfrequency omega. That one is
(1/2)*sqrt(4*c[0]*m[1]-d[0]^2)/m[1] for the ODE as written. And as d[0] is small it can be neglected so it is roughly sqrt(c[0]/m[1]) as you say.

Notice that even after doing  with(T); the result of dsolve(diff(x(t),t,t)+x(t)=0); is x(t) = _C1*sin(t)+_C2*cos(t), which means it is not expressed in terms of the T-versions but the global versions, which are available as :-sin and :-cos. Try
res:=dsolve(diff(x(t),t,t)+x(t)=0);
eval(res,t=360);
subs({:-sin = sin,:-cos=cos},res);
eval(%,t=360);

You may want to avoid redefining sin and cos and consider changing variables in the ODE, i.e. using another time tau:

PDEtools:-dchange({t=tau*k,x(t)=X(tau)},ODE,[tau,X]);
#And then pick k to fit your needs.
#If you use the angle in degrees as time tau, then you have tau = t*omega*180/Pi, so you should pick
k = Pi/(omega*180).
An additional comment. Why not just write s:=x; v:=D(x): a:=(D@D)(x):

If M > m >0 then it appears from the animation below that there is one real solution:
eq:=k = ln(( M-m)/m/(k-1));
res:=solve(eq,k);
#Animation ('eq' only depends on M/m, so we may safely use m=1):
plots:-animate(plot,[eval([k,rhs(eq),res],m=1),k=0..5],M=5..20);

If M < m there will be two real solutions when also M>0:

restart;
eq:=k = ln(( M-m)/m/(k-1));
res:=solve(eq,k,AllSolutions);
z:=op(indets(res,`local`));
res1:=subs(z=0,res);
res2:=subs(z=-1,res);
plots:-animate(plot,[eval([k,rhs(eq),res1,res2],m=1),k=-5..5],M=-5..10);
#This shows that there is tangency when k = 0 and M = 0:
solve({diff(eq,k),eq},{k,M});

It turns out that an explicit list is needed:

phi:=(x,y,z)->[2*x/(x^2+y^2+(-1+z)^2),2*y/(x^2+y^2+(-1+z)^2), 1+2*(-1+x)/(x^2+y^2+(-1+z)^2)];

In a procedure `transform/objects` local to the plottools package there is a conditional statement depending on the value of 'outdim' as defined by outdim := nops(f(i,j,outdim)); where f in our case is phi and the arguments to f are just (local) names. If phi(x,y,z) is seen as a sum of two terms, which it is in your original version, then outdim:=2.
However, if phi is defined as above, nops(phi(x,y,z)) returns 3 because there are 3 elements in the list.

In the 2-dim case there happens to be agreement between the two versions! Pure luck!

The example with Q below exhibits the problem. So you could define A implicitly or explicitly to be a table.

restart;
Q:=seq(a[k],k=1..3):
Q[1]:=3;
restart;
with(LinearAlgebra):
interface(rtablesize=infinity):
#In this loop A is implicitly defined as a table:
for k from 1 to 21 do A[k]:=Matrix(18) end do:
A[1]:=IdentityMatrix(18);
whattype(eval(A));

You should use 'add' instead of 'sum'.

'add' is for summations with concrete limits like in your case k = 0..3.

'sum' on the other hand tries to find a formula so that it actually might handle k=0..n, where n is not assigned a value. Also, and more importantly, 'add' and 'sum' have different evaluation rules: 'sum' evaluates its argument before acting, whereas 'add' does not. This means that MatrixPower(x,k) gets evaluated before k is known! That is the reason for the weird looking result you get. It is correct, though (at least in my example below).

Try the following:

with(LinearAlgebra);
A:=RandomMatrix(2);
c := 5;
x:=c*A;
MatrixPower(x,k);
'MatrixPower'(x,k);
sum(MatrixPower(x,k), k = 0..3);
simplify(%);
#You may or may not need evalf (A is randomly chosen)
add(MatrixPower(x,k), k = 0..3);
#Enclosing 'MatrixPower' in unevaluation quotes:
sum('MatrixPower'(x,k), k = 0..3);
#Notice that MatrixPower is more general than just using x^k (k concrete).
?MatrixPower
sum(x^k, k = 0..3);
add(x^k,k=0..3);
sum(MatrixPower(x,k),k=0..n);
eval(%,n=3);
simplify(%);
evalf(%);
#Compare with
add(x^k,k=0..3);

By enclosing the arguments in curly brackets you make them optional keyword parameters. In your first two calls to h1, there is no matching argument, so the default applies.

If you don't h1 to accept any other arguments than the ones explicitly given in the definition of h1, then end the argument sequence with a dollar sign:

h1 := proc( { [n,N] :: {positive,identical(n)=positive} := 1 }
         , { [l,L] :: {list,identical(l)=list} := [1,1,1] },$ )
  local %n,%l;
  if type(n,positive) then %n:=n; else %n:=rhs(n); end if;
  if type(l,list) then %l:=l; else %l:=rhs(l); end if;
      print([n,whattype(n),is(n,positive)]);
      print([%n,whattype(%n),is(%n,positive)]);
      print([l,whattype(l)]);
      print([%l,whattype(%l)]);
end proc :

I see nothing wrong with something like this (sd introduced only for a check).

Notice that the default doesn't have to have the declared type!

p:=proc({ [s,seed]::posint:= NULL} )
  local sd;
  if s<>NULL then sd:=randomize(s) end if;
  print([s,sd]);
  rand()
end proc;

p();
p(seed=45);

Does this do what you want?

restart;
g1:=proc({[n,N]::posint:=1})
  return g2(:-n=n);
end proc;
g2:=proc({n::posint:=1})
  return n^2;
end proc;

#And for your last question:

g3:=proc(n::{posint,identical(n)=posint}:=1)
local nn;
if type(n,posint) then nn:=n else nn:=rhs(n) end if;
return nn^2;
end proc;
g3(6);
g3(n=6);
g3();

If I understand your intention correctly, you want to turn the ode into a system of 3 first order odes.

If that is correct then you can use 'convertsys' from the DEtools package.

Due to an interesting bug (?) in DEtools, I begin by replacing the right hand side by just 'p(t)':

sys := diff(y(t),t$3) + 3*diff(y(t),t$2) + 5*diff(y(t),t)  = p(t);
#See the help page for convertsys for the arguments to convertsys. I chose xp as the name of the vector of derivatives.
DEtools[convertsys](sys,{},y,t,x,xp);
#Then replace 'p(t)' by the correct right hand side:
res:=subs(p(t)=diff(u(t),t$3) + diff(u(t),t$2) + diff(u(t),t),%);
#The system:
res[1];

#The translation (not exactly the one you mentioned, but probably what you meant):

res[2];

PS. I have submitted an SCR about the bug.

I suggest making the following changes besides ending the loop with a colon.

C:=Vector():
E:=Vector():

#Don't assign to i explicitly, but let the loop do it (and start from 1)

for i from 1 while t.....

C(i):=t:  #Notice the use of a parenthesis, not square brackets: Programmer indexing
E(i):=p_2:  #Same thing
#Again: Don't assign to i. The loop does it.

Outside the loop:

with(plots):
pointplot([C,E]);

#You could also (if you like) tighten up the definitions of S, A_zu, and A_ab a little (inside the loop):

S:=piecewise(
   t <= 9.5,s_1/1000,
   t <= 12.875,s_2/1000,
   t <= 15.5,s_3/1000,
   t <= 20.75,s_4/1000,
   t <= 26.75,s_5/1000
                       );
if S >= 0.01 then A_zu := 0 else A_zu := b*h-b*S end if;
if S < 0.01 or S >= 0.018 then A_ab := b*h-b*S else A_ab:=0 end if;

It appears that you have a couple of errors in your procedure: (1) a missing 'then' , (2) a missing condition in 'if lambda' (presumably you mean 'if lambda < lambda0'), and (3) 'e^(...)' should be 'exp(...)'.
I had problems in copying my corrected procedure into MaplePrimes. I suspect (1) and (2) are actually MaplePrimes copying mistakes.

There is no need for Cp (it doesn't hurt, though).

So maybe this is what you intended:

Cplot := proc (beta, lambda)
local lambda0;
lambda0 := 0.22e-1*beta^2+(7/1000000000000000)*beta+5.6;
if lambda < lambda0 then
   0.1e-2
else
  ((1/2)*lambda+(1/2)*(-1)*0.22e-1*beta^2+(1/2)*(-5.6))*exp((-1)*.17*lambda)
end if
end proc;

#Testing:
Cplot(5,40);
#But symbolic input gives an error:
Cplot(a,b);
#Use unevaluation quotes to prevent that
'Cplot'(a,b);
#So here is the plot:
plot([seq('Cplot'(k*5,lambda),k=0..6)],lambda=0..60);

#An alternative to using unevaluation quotes outside the procedure is to take care of the problem inside as is done in this version:
Cplot := proc (beta, lambda)
local lambda0;
if not type([beta,lambda],list(realcons)) then return 'procname(_passed)' end if;
lambda0 := 0.22e-1*beta^2+(7/1000000000000000)*beta+5.6;
if lambda  < lambda0 then
   0.1e-2;
else
   ((1/2)*lambda+(1/2)*(-1)*0.22e-1*beta^2+(1/2)*(-5.6))*exp((-1)*.17*lambda)
end if
end proc;

plot([seq(Cplot(k*5,lambda),k=0..6)],lambda=0..60);

We have no way of actually checking your worksheets since we don't have the matrices

PsiST := ImportMatrix("L:/MapleTabellen/[PsiTST]_NT=5e17,xT=100,ET=1700.dat", source = delimited);

FD32 := ImportMatrix("L:/FD/FD32_M.txt", source = delimited);

From a very superficial look at your worksheets it appears that the value of Digits is the default, i.e. 10.

You may want to raise that.

If you use square brackets around the N[i]'s then the order is respected:

solve(List,[seq(N[i],i=1..9)]);
op(op(%));
assign(%);
List;

Warnings are printed with 'print', see

showstat(WARNING);

so are not output.

Here are two ways of saving the warning messages for inspection or handling later.

The first one redefines the WARNING procedure to make it assign to the global variable 'lastwarning'.
The second method makes use of a global variable `DSI/Warnings`, which is used by 'odeplot' when _Env_in_maplet is true and is assigned the list of warnings issued. I don't use maplets so I have no idea if the assignment _Env_in_maplet:=true will have undesirable consequences in your case.

restart;
showstat(WARNING);
unprotect(WARNING);
WARNING:= proc(msg::{string, symbol}) global lastwarning;
   lastwarning:=msg;
   print(INTERFACE_WARN(1,msg,args[2 .. -1]))
end proc;
protect(WARNING);
sol:=dsolve({diff(y(x),x)=y(x)^2,y(0)=1},numeric);
plots:-odeplot(sol,[x,y(x)],-1..1);
lastwarning;
lastexception;
sol(1);
lastexception;
StringTools:-FormatMessage(lastexception[2..-1]);
WARNING(%);
lastwarning;
#Instead of redefining 'lastexception' at each call to WARNING, 'lastexception' could be assigned a cumulative list of all warnings just as is done by `dsolve/numeric/warning` in the next method: 
restart;
#showstat(`plots/odeplot`);
showstat(`dsolve/numeric/warning`);
_Env_in_maplet:=true;
sol:=dsolve({diff(y(x),x)=y(x)^2,y(0)=1},numeric);
plots:-odeplot(sol,[x,y(x)],-1..1);
`DSI/Warnings`;