Maple 2018 has the ability to encrypt procedures. In "What is new"  it says:
"Maple 2018 lets you encrypt procedures. An encrypted procedure acts and behaves like any normal procedure in Maple, but it can not be viewed or debugged. This provides a way for professional users to share executable procedures while protecting their IP."

The result returned by int with fc is wrong.
 

restart;
f := ((1/2-I*t)^(-s)-(1/2+I*t)^(-s))/(2*I);
fc:=evalc(f);
res1:=int(f,t=0..infinity) assuming s>1;
int(fc,t=0..infinity) assuming s>1;
res2:=simplify(%);
# tests:
eval([res1,res2],s=2); # result [2,8]
evalf(Int(eval(fc,s=2),t=0..infinity)); # 2
evalf(Int(eval(f,s=2),t=0..infinity)); # 2
## Seeing all results:
int(fc,t=0..infinity,method=_RETURNVERBOSE) assuming s>1;
int(f,t=0..infinity,method=_RETURNVERBOSE) assuming s>1;

So for fc the "successful" methods are ftoc and ftocms and they agree and are wrong.
For f the "successful" methods are ftoc and ftocms as for fc, but in addition the correct meijerg.
## NOTE added: Today I cannot reproduce the wrong ftoc and ftocms results for f itself. They now  agree with the correct  meijerg result.

The simple exceptional solution not found by dsolve(ode) is found by `dsolve/IC/easy` when dsolve is given the initial value problem.
It can be found from the generic initial value problem by taking a limit:
 

restart;
ode:=diff(y(x),x)=x*ln(y(x)):
ic:=y(1)=y0: # generic ic
sol:=dsolve({ode,ic});
map(limit,sol,y0=1); # OK
## To see that the original initial value problem is handled by `dsolve/IC/easy`:
debug(`dsolve/IC/easy`);
dsolve({ode,y(0)=1});

Occasionally simplify makes an expression less simple than it was before or just not as simple as you had hoped.
In situations like that it can be worth while to try using map(simplify, expr).
Examples:
 

restart;
A := (1 - cos(s)^2)/sin(s)^2;
B := (1 - cos(t)^2)/sin(t)^2;
map(simplify,A*B);
##
S:=2*A+4*B;
simplify(S); # not simple
map(simplify,S); # 6

Like tables, procedures, and modules, records have last name evaluation.
Example: sin; evaluates to just sin, i.e. to the name of the sin procedure.
Try:
 

r:=Record(a,b);
type(r,record);
whattype(r);
whattype(eval(r));
whattype(sin);
whattype(eval(sin));

There is a help page:
? last name evaluation
#################
About the acceptance of strings:
Not accepted in Maple 15 (and earlier at least not in Maple 12).
But accepted from Maple 16.02 and up it seems.
In Maple 16.02 and also in the present Maple 2018.1 we find in the help page for Record the statement:

"To make it easier to construct records programmatically, the symbol in a record field specifier may also be passed to the Record constructor as a string. This avoids potential evaluation of field names that correspond to an assigned global variable."

You don't give the details about the solve statements so I shall just give you a rather trivial example.
(Obviously you are missing an 'and' in your code.)
 

restart;
bound1:=NULL; bound2:=45;
if bound1<>NULL and bound2<>NULL then blahblah else booh end if;
bound1:=NULL; bound2:=NULL;
if bound1<>NULL and bound2<>NULL then blahblah else booh end if;
bound1:=78; bound2:=45;
if bound1<>NULL and bound2<>NULL then blahblah else booh end if;

Using method = laplace for linear systems with constant coefficients often leads to simpler looking results.
I used the code for the system typed by mmcdara.
 

restart;
sys := {
diff(tau__1(t),t)=phi__1(t),
diff(tau__2(t),t)=phi__2(t),
diff(phi__1(t),t)+2*phi__1(t)+3*phi__2(t)+4*tau__1(t)+5*tau__2(t)=6,
diff(phi__2(t),t)+8*phi__1(t)+7*phi__2(t)+9*tau__1(t)+10*tau__2(t)=11
};
## The dependent variables have to be given explicitly when using method=laplace:
fus:=indets(sys,And(function,Not(specfunc(diff))));
res:=dsolve(sys,fus,method=laplace);
fus0:=subs(t=0,fus);
ics:=fus0=~{1,-1,2,-2}; #Example
res2:=allvalues(res):
plot(subs(res2,ics,fus),t=0..2);
## A floating point version:
res3:=simplify(fnormal(evalf[15](res2),9),zero); # floating point version
plot(subs(res3,ics,fus),t=0..2);

The first plot:

Note. You may want a legend. Thus you need a list of functions, not a set.
Modify the plot command as follows:
 

fuslist:=convert(fus,list);
plot(subs(res2,ics,fuslist),t=0..2,legend=fuslist);

Simplifying ode helps:
 

restart;
F := x * ( y(x) + x*sqrt(x*y(x)) + sqrt(x^3*y(x)) );
ode:= diff(y(x),x) = F;
ode1:=simplify(ode) assuming positive;
dsolve(ode1);

Continuing with an explicit solution:
 

dsolve(ode1,explicit) assuming x>0;
allvalues(%);

############################
These lines may explain what causes the error (?):
 

infolevel[dsolve]:=5:
ode2:=PDEtools:-dchange({y(x)=v(t),x=sqrt(t)},ode,[v,t]);
ode3:=simplify(%) assuming t>0 ;
dsolve(ode2); # A different error, but still ...
dsolve(ode3); # OK

Here is a version that uses dsolve/numeric.
 

restart;
Omega := 2*Pi*N; 
R0 := a*tanh((a^2-mu)/(2*T_c))*ln((2*a^2+2*a*q+q^2-2*mu-I*Omega)/(2*a^2-2*a*q+q^2-2*mu-I*Omega))/q-2;                      
T_c := 0.169064e-1; 
mu := .869262; 
#N := 10;
#R1 := int(R0, a = 0.1e-3 .. 100);    
#R2 := evalf(abs(R1));

ode:=diff(A(a),a)=-R0;
res:=dsolve({ode,A(100)=0},numeric,parameters=[N,q],output=listprocedure);
R:=subs(res,A(a));
Q:=proc(N,q) 
  if not type([N,q],list(numeric)) then return 'procname(_passed)' end if;
  res(parameters=[N,q]); 
  abs(R(1e-4)) 
end proc;
Q(10,5); #Test
t0:=time():
plot(curry(Q,10),1e-3..10,caption="N = 10");
time()-t0; # 28s on my machine
t0:=time():
plots:-animate(plot,[Q(N,q),q=1e-3..10],N=0..10,frames=11); #Only 11 frames
time()-t0; # 350s on my machine

I should point out that acer's method is much faster, also on my machine. One plot by his method takes about a second, where the one plot above takes about 30 seconds.

Maybe you want the coefficients in a taylor expansion in p?
I'm using mtaylor instead of taylor for convenience since it omits the big O term:
So you can do:
 

T_HPMEq1:=mtaylor(HPMEq1,p,N+1):
T_HPMEq2:=mtaylor(HPMEq2,p,N+1):
T_HPMEq3:=mtaylor(HPMEq3,p,N+1):
#
for i from 0 to N do
equ1[i] := coeff(T_HPMEq1, p, i) = 0; 
equ2[i] := coeff(T_HPMEq2, p, i) = 0;
equ3[i] := coeff(T_HPMEq3, p, i) = 0 
end do;

Since odetest handles the implicit given solution nicely and regardless of the name or value of the arbitrary constant (see my reply to Carl Love) another workaround is to convert the "explicit" solution to an implicit one. I use quotes since a solution having a RootOf can hardly be called a genuinely explicit one.
 

restart;
ode := diff(y(x), x) = x*ln(y(x)):
sol:=dsolve(ode, implicit);
odetest(sol,ode); #OK
odetest(subs(_C1=C1,sol),ode);#OK
odetest(subs(_C1=1,sol),ode); #OK
explicit_sol:=op(solve(sol,{y(x)}));
DEtools[remove_RootOf](subs(_C1=C1,explicit_sol));
odetest(%,ode);
DEtools[remove_RootOf](subs(_C1=1,explicit_sol));
odetest(%,ode);

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]);