For your example

fileName := "C:\\foo.txt";
fd       := fopen(fileName,WRITE):

There are two possible outcomes, either

1. fopen() is successful so carry on. Why/what would you check? or
2. fopen() fails and Maple generates an Error

As a rule of thumb, whenever any software detects an error, it stops with an (un)informative error message. If you plan on continuing afer an Error has been generated, then it is your responsibility to "trap" the error and decide what to do next. Consider the foloowing example of opening a non-existent file in READ mode.

   restart;
   fname:="D:/Users/TomLeslie/myMaple/test1.txt":
#
# Verify that above filename does not exist
#
   FileTools:-Exists(fname);
   try
       #
       # So if one "opens" a non-existent file
       # for read, then an ERROR is to be
       # expected
       #
          f:=fopen(fname, READ)
  catch:
          #
          # Catch any and all errors. Print a polite warning
          #
             printf("Ooooops! something went wrong with fopen\n");
          #
          # So add more code here to do something different
         #
    end try;

which produces
                             false   <- from the FileTools:-Exists() statement
                            Ooooops something went wrong with fopen <- from the catch clause

You could do the whole thing analytically - no need to use dsolve(...numeric).

Makes coding simpler, runs much faster, it is much easier to manipulate the solutions of the ODEs, and it seems to come up with the same answer as the 'numerical' version. See attached

ODEplot3.mw

I managed to get a solution for your second case using the DirectSearch() addon package, available from the Maple Application centre.

See the final group of the attached worksheet for the answer I obtained. I have saved/upoaded this worksheet with output, since I am assuming that you will not be able to re-execute it (because you probably don't have the DirectSearch package installed)

bigSys.mw

From the information you supply, one can plot any function of quantities returned by the dsolve() commands.

For illustration, the attached plots x[i](t) versus t, for i=1..5

ODEplot.mw

In the code you supplied, there is a 'curvefitting' execution group whose output doesn't seem to be used for anything. Also Phi is undefined, and I have no idea whether this is important or not

Acer and vv have already answered this problem - but just because alternatives are good (and stitching everything together)

  restart;
  with(LinearAlgebra):
  k__1 := (12*2)*10^6/3^3:
  k__2 := k__1:
  k__3 := (3*1.5)*10^6/3^3:
  k__4 := (12*1.5)*10^6/3^3:
  k__5 := 12*10^6/3^3:
  m__1 := 8000:
  m__2 := 7000:
  m__3 := 6000:
  m__4 := 5000:
  m__5 := 10000:
  K := Matrix(5, 5, [ [  k__1+k__2,  -k__2,             0,                    0,                   0       ],
                                [ -k__2,             k__2+k__3,  -k__3,              0,                   0       ],
                                [  0,                  -k__3,             k__3+k__4,   -k__4,             0       ],
                                [  0,                   0,                  -k__4,              k__4+k__5,  -k__5 ],
                                [  0,                   0,                   0,                   -k__5,             k__5 ]
                             ]
                     ):
   M := Matrix(5, 5, [ [m__1, 0,       0,          0,           0          ],
                                  [0,       m__2, 0,         0,            0         ],
                                  [0,       0,        m__3,  0,            0         ],
                                  [0,      0,         0,         m__4,     0         ],
                                  [0,      0,         0,         0,            m__5  ] 
                               ]
                    ):
   E,Q:=map(simplify, [Eigenvectors(K, M)], zero)[];
   perm:=sort(E,output=permutation):
   E:=E[perm];
   Q:=Matrix([seq(Q[..,perm[j]], j=1..5)]);

Would that be y(t), diff(y(t),t), or some other quantity????

The following will return the maximum values of y(t) and diff(y(t),t) over the specified range

  restart;
  res := dsolve(  { 25*(diff(y(t), t, t))+4*(diff(y(t), t))-3*y(t) = cos(3*t),
                              y(0) = 0,
                              D(y)(0) = 1
                            },
                           numeric,
                           output=listprocedure
                        );
  plots:-odeplot( res,
                             [  [t, y(t)],
                                [t, diff(y(t),t)]
                            ],
                            t=0..4
                         );
  with(Optimization):
  Maximize(rhs(res[2]), 0..4); # <- return maximum of y(t)
  Maximize(rhs(res[3]), 0..4); # <- return maximum of diff(y(t),t)

but you can 'isolate' anyy of the variables in your expression with

eq:=(N*(P-p__th)/K+p__th)/(2+2*N*T/K+(1+N*T/K)^2/(N*(P-p__th)/K+p__th))=p__th/(2+1/p__th);
seq( isolate(eq, kk), kk in indets(eq));

setting Digits=50 at the start of your worksheet - works for me

You probably ought to consider why a high numerical accuracy is often required when dealing with the numerical integration of a an oscillatory function

1. The worksheet you supplied was doing a few "odd" things when I started to debug. Not necessarily wrong - just odd. So I ened up converting/rewriting a lot in 1-D math input. Ended up writing everything in 1-D math
2. To assist debug I ended up defining all quatities which contained the independent variable ('y' in your case) as 'functions' of this variable. This made it easier to check whether or not these functions were reasonably well-behaved
3.  At the same time I converted all of jour "indexed" variables, to 'subscripted' variables, rg converted 'w[s]' to 'w__s', Partly becuase you are never using these as 'indexed' variables and partly because the evaluation rules are slightly different for 'indexed' variable uses in function definitions, and I wanted to avoid, loads of extra 'eval()' statements
4. So far, nothiing I have done should have changed the "logic" of your calculation. All of the arguments to the ellitpicPi functions appear to be reasonably well-behaved - except maybe when the independent variable is <~0.1, where some complex values occurred. However in your final definition of 'alpha' you use a mixture of 2-argument' and 3-argument ellipticPi functions. The original 3-argument ellipticPi() functions return complex values for nearly all values of the independent variable, except for a small region around 0.1. However when I examine these 3-argument ellipticPi functions, the order of the arguments seems a little suspicious, when compared with 2-argument ellipticPi functions and the Maple help page. With no justification at all, I flipped the order of the first two arguments in the 3-argument ellipticPi() functions. Only you know whether this was a typo in your original
5.  With all of the above, and assuming I haven't  made a typo somewhere, the attached at least produces an answer

elipInt.mw

works for me

  F:= t-> piecewise( t < 0, 0,
                                 t < 1, 100,
                                 t < 2, 200,
                                 t < 3, 50
                              );
  sol:= dsolve( [  diff(y(t),t$2)+4*diff(y(t),t)-3=F(t),
                             y(0)=0,
                             D(y)(0)=1
                         ]
                      );
  plot(rhs(sol), t=0..5);

Use the following code, with 'fname' set to appropriate something for your machine

Download ArrXLexp.mw

I haven't been able to make option  3 work - I'm guessing that this i because Excel allows multiple entries in a given cell, these have to be "structured" in  some specific way in an xlsx file - and I can't figure this out

Since x(t) and diff(x(t), $n) are both of type  'function' and both of type 'dependent', I'm not sure that you can distinguish these with 'indets()' - although I might be wrong. Furthermore, the independent variable is of type 'name'.

The following seems to work, but is really 'ugly'

f:= x-> `if`( StringTools:-Search("diff(", x)=0,
                   parse(x),
                   NULL
                ):
g:= x-> map( convert,
                      `union`
                      ( indets(eq, name),
                        indets(eq, function)
                      ),
                     string
                   ):
#
eq:= diff(x(t), t) = x(t) + 1  ;
f~(g(eq));

It is well worth reading the help at ?last name evaluation - it is occasionally subtle.

As has been mentioned growing lists in this way is a bad habit, since lists are not mutable. For small examples, it really doesn't matter much, but it is a bad habit to get into - I know, beauase I used to do this pretty regularly,  and eventually I got "bitten". For your 'toy' example a better way would probably be

restart;
nonIdMaps := [ seq
                           ( `if`
                             ( f(y) <> y,
                               eval(f),
                               NULL
                           ),
                           f in [ x -> x,x -> 2*x,x -> 3*x ]
                         )
                      ];

Use the menu sequence

Format-->Styles->Maple Plot->Modify->Alignment

and toggle this to "Left" rather than "Center"

The following will provide *an*  answer in Maple2017.1

sol:=dsolve( { (x-1)*diff(y(x), x$2)+diff(y(x),x) = 1+x,
                         y(0)=0,
                         y(1)=1},
                         type = numeric,
                         method = bvp[midrich],
                        abserr=10000
                     );
plots:-odeplot(sol, [x,y(x)], x=0..1);

The "discontinuoty" at x=1 is pretty obvious. It would seem that this boundary condition is not really compatible with the condition y(0) and the ODE. You may want to reconsider