Use cat.
 

x1, x2, x3, x4:=1,2,3,4;
for i from 1 to 4 do unassign(cat(x,i)) end do;
x1,x2,x3,x4;

As for || see the help page ?|| where it says:
"Note: The use of this operator is discouraged. Where possible the function cat should be used."

I shall just address the warning about _EnvDSNumericSaveDigits. That is curious and I haven't seen that before.
It seems to come from `dsolve/numeric` or rather from one ot its auxiliary procedures.
Notice that the warning disappears if you make an explicit integer assignment to _EnvDSNumericSaveDigits such as
 

_EnvDSNumericSaveDigits:=15:

Since the environmental variable _EnvDSNumericSaveDigits is meant to save the old value of Digits any positive integer will make the warning go away.
The warning must come up in in your worksheet because somehow _EnvDSNumericSaveDigits avoided getting assigned to in your worksheet. The warning doesn't come up by a simple call to any of the two procedures, only in the fsolve context (or so it seems).
### Furthermore if you remove the range argument then the warning disappears.
Here is a much simpler example where the same warning comes up:
 

restart;
ode := diff(x(t), t) = x(t);
t1:=2:
h := proc (x0)  
  local X;    
  X := eval(x(t), dsolve({ode, x(0) = x0}, numeric,range=0..t1, output = listprocedure)); 
  X(t1)-1
end proc:
h(0);
fsolve(h);
h(%);

As you see the result is correct, but the warning is confusing.
Try removing the range argument and the warning is gone.
I would have done this example using the parameters option in this way:
 

restart;
ode := diff(x(t), t) = x(t);
t1:=2:
res:=dsolve({ode, x(0) = x0}, numeric,range=0..t1, output = listprocedure,parameters=[x0]);
X:=subs(res,x(t));
h := proc (x0)  
  res(parameters=[x0]);   
  X(t1)-1
end proc:
h(0);
fsolve(h);
h(%);

No warnings. You can remove the range argument (and still no warnings).
I shall submit an SCR (bug report) about the warning as it is confusing to the user.

If you save both as .mpl files and use the complete path and then use a read statement to read main_module.mpl your setup with $include should work.
I just tried:
 

main_module:=module()
 option package;
 local C; 
 local sub_module;
 export main_entry; 

 C:=99; #see if this can be "seen" from child module
	
$include "F:/MapleDiverse/MaplePrimes/sub_module.mpl"

 main_entry :=proc()
   print("in main_module:-main_entry()"):
   sub_module:-main_entry();
 end proc;
end module;

and
 

#sub_module.mpl
sub_module:= module()
    export main_entry;
    main_entry :=proc()
     print("In main_module:-sub_module:-main_entry(), C=",C):
    end proc;
 end module;

The files are in my case both in the folder "F:/MapleDiverse/MaplePrimes".
Now in Maple use the read statement
 

restart;
read "F:/MapleDiverse/MaplePrimes/main_module.mpl";
## Test:
main_module:-main_entry();

It appears that Maple 2018 can find a symbolic solution. It is big, but integrating the first ode w.r.t. x and using ics helps and is done here:
Note added: The code below assumes that Q is a constant. If Q depends on x the method may still work. If not use numerical solution as Carl does.

restart;
sys_ode:=    2*C__5*diff(w(x),x) - C__1*diff(u(x),x$2) - 2*C__4*diff(w(x),x$3) = Q, 
   2*C__2*u(x) - C__1*diff(w(x),x$3) - 2*C__3*diff(u(x),x$2) = 0;
ics:= w(0) = 0, D(u)(0) = 1, (D@@2)(w)(0) = 0, D(u)(100) = 1, (D@@2)(w)(100) = 2;
###
map(int,sys_ode[1],x=0..x,continuous);
ode1:=eval(%,{ics});
ode2:=sys_ode[2];
### Now solve symbolically:
res:=dsolve({ode1,ode2,ics}):
U,W:=op(subs(res,[u(x),w(x)])):
Digits:=200:
## Example plot:
plot(eval(U,{Q=1,C__1=1.2e8, C__2=2e9,C__3=2e8,C__4=0,C__5=3e7}),x=0..100,size=[2000,default]);

The plot of W is done similarly:
 

plot(eval(W,{Q=1,C__1=1.2e8, C__2=2e9,C__3=2e8,C__4=0,C__5=3e7}),x=0..100,size=[2000,default]);

In procB you have the formal parameter the_equation. That appears in the procedure body of procB in the print statement as well as in the call to procA. Wen procB is called with the equation y=3 the formal parameter the_equation will be replaced by y=3 in all its occurrences in the body of procB. Thus the call to procA will be procA('y=3' = (y=3)), so that the default NULL will be used by procA.

You can use this instead:
 

procA(':-the_equation'=the_equation);

Your sys1 will necessarily have a singularity for some t1 < tini.
The reason is that as long as a(t) > 0 it follows that a(t) is increasing. Going back in time a(t) will decrease to 0 at some finite value of t (called t1 above). Here a'(t1) = infinity.
To prove that that is what happens is rather straightforward. Your ode is of the form
 

ode:=diff(a(t),t) = k*a(t)*sqrt(1/a(t)^3 + m); # k>0, m>0.

with a(t0) = a0 > 0.

As long as a(t) > 0 we have
 

diff(a(t),t) > k*a(t)*sqrt(1/a(t)^3) = k/sqrt(a(t))

Rewrite this to get
 

sqrt(a(t))*diff(a(t),t) > k

Integrate from t to t0 (with t < t0) to get
 

2/3*(a0^(3/2) - a(t)^(3/2)) > k*(t0-t)

from which follows that
 

t > t0 - a0^(3/2)/k*(2/3)

since we are assuming that a(t) > 0. Thus a(t) will necessarily approach 0 from above at a finite value of t < t0.

Try this where I have replaced the letter i by the letter I:
 

restart;
A:=Matrix(6,(i,j)->cos(2*Pi/6*(i-1)*(j-1))-I*sin(2*Pi/6*(i-1)*(j-1)));
B:=conjugate~(A);
expand~(A.B);

If you eliminate one of CNO and CNH3 you can easily solve the system. From your two odes it follows that the difference between CNO and CNH3 is constant, thus in the case presented equal to zero.
 

restart;
r := 2900*exp(-53300/(R*T))*CNO(t)^0.62*CNH3(t)^(-0.05);
ode := diff(CNO(t), t) = -r: #(negative because they are decaying with time)  
ode2 := diff(CNH3(t), t) = -r:
ics := CNO(1020) = 1.6, CNH3(1020) = 1.6; #(sets up known initial conditions)
ode := subs(R = 8.314, T = 473, ode):
ode2 := subs(R = 8.314, T = 473, ode2):
#sys_ode := (ode, ode2) ;
odeS:=subs(CNH3(t)=CNO(t),ode);
#dsolve({odeS,CNO(1020)=1.6}); #??? Why does it take so long (or forever?)
## So do instead:
sol:=dsolve(odeS);
eval[recurse](sol,{t=1020,CNO(1020)=1.6});
solve(%,{_C1});
eval(sol,%);
res:=solve(%,CNO(t)); # The solution for CNO(t)
t0:=solve(res=0,t);
plot(res,t=1020..t0);

The solution becomes 0 in finite time t0 (approximately 1775.4) due to the power of the right hand side of odeS being positive and less than 1. After that the solution continues as zero.

If the initial conditions are different it appears that you run in to an integral that Maple cannot find. A numerical solution may then be used:
 

ics2 := CNO(1020) = 1.6, CNH3(1020) = 1; #Example of different initial conditions
# Now CNH3(t)-CN0(t) = constant = 1 - 1.6 = -0.6
odeS2:=subs(CNH3(t)=CNO(t)-0.6,ode);
res:=dsolve({odeS,CNO(1020)=1.6},numeric);
plots:-odeplot(res,[t,CNO(t)],1020..2000);

Notice though that since CNH3(t) = CNO(t) - 0.6 we must stop the integration when CNH3(t) = 0.
You could just plot CNH3(t) like this:

plots:-odeplot(res,[t,CNO(t)-0.6],1020..2000);

An alternative:

resE:=dsolve({odeS,CNO(1020)=1.6},numeric,events=[[CNO(t)-0.6,halt]]);
resE(2000);

You will be told that integration stopped at t = 1279.9381.
 

Make the following addition and change:
 

col:=["black", "blue", "red", "pink"];
for k to 4 do 
   R := dsolve(eval({bc, sys}, Ha = L[k]), fcns, type = numeric, AP); 
   AP := approxsoln = R; 
   p1u[k] := odeplot(R, [y, U(y)], 0 .. 1, numpoints = 100, labels = ["y", "U"], style = line, color = col[k]) 
end do:

Begin by revising the definition of F and then use a revised version of your second method.
## I revised my answer too  :-)
## Introduce Y'' as Yd2(x) (the easiest solution in this case, I think):

F := dsolve({cond, sys,diff(Y(x),x,x)=Yd2(x)}, numeric, output = listprocedure);
Y1 := subs(F,Y(x));
Y2 := subs(F,diff(Y(x), x));
Y3 := subs(F,Yd2(x));

The remaining ones ( R(x), mu(x) , and mu'(x) ) are found quite similar to Y(x), and Y'(x).

## There appears to be singularity problems.
 

F := dsolve({cond, sys}, numeric, output = listprocedure); #No Yd2
plots:-odeplot(F,[x,Y(x)],0..Pi);

You seem to be interested in x = 0..Pi, but you get from odeplot:

Warning, cannot evaluate the solution further left of 2.6226594, probably a singularity

Furthermore when trying to include diff(Y(x),x,x) = Yd2(x):
 

F := dsolve({cond, sys,diff(Y(x),x,x)=Yd2(x)}, numeric, output = listprocedure);
plots:-odeplot(F,[x,Y(x)],0..Pi);

then we get an error I don't remember having ever seen before:

Error, (in dsolve/numeric/SC/IVPrun) invalid reversal of integration direction

Going the other way is OK at least until 3.66:

plots:-odeplot(F,[x,Y(x)],Pi..3.66);

# When I use the second version of F (i.e. the one including Yd2) and try F(3) then I get this error:

Error, (in unknown) cannot evaluate the solution further left of 3.6605260, probably a singularity

which is surprising to me because who told F to go from Pi to 3.6605260 before going to 3?

## I notice that the weird error disappears if cot(x) is replaced by 1. Singularity problems don't disappear though.
###############
The problem (unrelated to other singularities) appears to have to do with the initial singularity at x = Pi where cot(x) is infinite (undefined):
Try

F := dsolve({cond, sys}, numeric, output = listprocedure,abserr=1e-15,relerr=1e-13);
F(3);

and you get the error message:

Error, (in unknown) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

Actually, in Maple 2018 I get the complaint:

Warning, cannot evaluate the solution further right of 3.4774930, maxfun limit exceeded (see ?dsolve,maxfun for details)
So since maxfun can be set, I did this:
 

res := dsolve(sys union ics, numeric, maxfun=10^6);

But the complaint just changed the t-value slightly to 3.4775188.
There may be an actual singularity though.
Try this:
 

sys:=solve(odes,diff~({x,y,z}(t),t)); # Solve for the derivatives
## Then do
plots:-odeplot(res, [t,rhs(sys[1])], Pi .. 4,caption="The derivative of x(t)");

What I see is this:

Using method=rosenbrock gives me the error you report:
 

res := dsolve(sys union ics, numeric, maxfun=10^6,method=rosenbrock);

The correct way to define a function is like this:
 

T:=x->1/2*x*sqrt(-x^2+(60-x)^2);
##
solve(diff(T(x),x)=0,x);
##
T(20);

If you haven't already done so, you should take a look at Student:-Calculus1:-ApproximateInt and Student:-NumericalAnalysis:-Quadrature at least for comparison with your own code.
 

restart;
fx:=x^2*sqrt(25-x^2);
Student:-Calculus1:-ApproximateInt(fx, x=0..5,method=trapezoid,partition=100);
evalf(%);
Student:-Calculus1:-ApproximateInt(fx, x=0..5,method=midpoint,partition=100);
evalf(%);
Student:-NumericalAnalysis:-Quadrature(fx,x=0..5,method=trapezoid, partition=100);

## Note: The following NumericalAnalysis:-Quadrature method is the same as the Calculus1:-ApproximateInt midpoint method:
 

Student:-NumericalAnalysis:-Quadrature(fx,x=0..5,method=newtoncotes[open,0], partition=100);

mmcdara already pointed out that you are missing the integration variable in the procedure pn.
But pn has another problem: pn has the formal parameter t  ( pn:=proc(i,n,t) ... ). OK.
But in the body of pn you have int(hd[i],t). That forces you to use literally t when using the procedure because hd[i] depends on the global t. This is bad (weird):
 

hd:=t^2:
p:=proc(t) 
int(hd,t)
end proc:
p(t);
p(s);

Simple solution for this case: leave out the third formal argument in pn.
##
You write:  " #alpha1(t) and alpha2(t) must be defined as functions of t ".
but you don't do that.
And what is f?
(I agree with mmcdara that you should use add instead of sum for concrete finite sums).
###########################
Putting in the functions from Example 2 in your reference just before defining RHS:
 

f:=t->t^2;
alpha1:=t->-1:
alpha2:=t->t:
y0:=1:
y1:=-1:

Then run your corrected worksheet.
After that do for comparison with the numerical solution found by Maple (which is good):
 

ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2;
res:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric);
plots:-odeplot(res,[t,Y(t)-y(t)],0..1);

The error in the method with the given choices:

With J=7 (so N=128) you get the error:

When doing this:
 

restart;
A:=LinearAlgebra:-RandomMatrix(5);
B:=LinearAlgebra:-RandomMatrix(5,datatype=float);
C:=LinearAlgebra:-RandomVector(50,datatype=float);

I see in the Variables list to the left:
A   5x5 Matrix
B   5X5 Matrix
C   50-element Vector[column]

###
Adding an Array:
 

E:=Array(1..100,i->i);

makes me see also
E   100-element Array.
###
`Standard Worksheet Interface, Maple 2018.0, Windows 10, March 10 2018 Build ID 1298750`