When using interface(typesetting=standard, the default in Maple 2015) the printing is done without parentheses.
In Maple 2018 the default setting of typesetting is extended.
So in Maple 2018 and Maple 2015 try:
 

restart;
interface(typesetting=extended); # The default in Maple 2018
t := table([1=table([a=123]) ]);
save t,"F:/MapleDiverse/test.m";
####
restart;
interface(typesetting=standard); # The default in Maple 2015
read "F:/MapleDiverse/test.m";
t[1][a] := 321;

In the first part (the saving) you can also try using interface(typesetting=standard) instead. It doesn't make any difference for the second part as long as that is done with typesetting=standard.

So my conclusion is that it just is a typesetting matter.

There has been a change from Maple 2015 to Maple 2016.
I tried the following in Maple 12, 17, 2015, 2016, 2017, and 2018:
 

restart;
type(fsolve,`module`);

In Maple 2016 and later fsolve is a module having ModuleApply as a local. Thus fsolve can be used as if it were a procedure.
In earlier version fsolve was just a procedure.
Incidentally, didn't you mean reinitialize instead of initialize?
If you try forget(fsolve,xxx=true) you don't get a complaint, thus a wrong optional argument of the type name= anything is just ignored.
The help page for forget has this statement about reinitialize:

"reinitialize = true or false (default: true)
This option only applies to the case when f is a procedure which is not a member of a module. In this case, if reinitialize=true, the corresponding procedure f is completely discarded and reloaded from its Maple archive, causing its remember table to be restored to its initial state." (my italics).

Try this:
 

op(2,eval(forget)); # forget is also a module
showstat(forget::ModuleApply);

Remove the range argument to dsolve and then do SOL(-2);
That keeps you from getting confused.
With the range argument just do SOL(-2);  before you go on.

As far as I can see your boundary value problem has no solution.
 

restart;
L:=75/10^6; # rational
eq1:=0.2079268293e-3*(diff(psi(x), x, x))-101.2059621*(diff(u(x), x, x))-101.2059621*(diff(w(x), x))*(diff(w(x), x, x));
##
eq2:=0.2079248162e-3+29.57317072*psi(x)-29.57317072*(diff(w(x), x))-1.996189024*10^(-9)*(diff(psi(x), x, x))-1.996189024*10^(-9)*(diff(w(x), x, x, x));
##
eq3:=151.8089432*(diff(w(x), x))^2*(diff(w(x), x, x))+101.2059621*(diff(u(x), x, x))*(diff(w(x), x))+101.2059621*(diff(u(x), x))*(diff(w(x), x, x))-0.2079268293e-3*(diff(psi(x), x, x))*(diff(w(x), x))-0.2079268293e-3*(diff(w(x), x, x))*(diff(psi(x), x))+29.57317070*(diff(w(x), x, x))-29.57317072*(diff(psi(x), x))-1.996189024*10^(-9)*(diff(psi(x), x, x, x))-1.996189024*10^(-9)*(diff(w(x), x, x, x, x))-p*(diff(w(x), x, x));
sys:=convert({eq1,eq2,eq3},rational): # for casesplit
bcs:=psi(0) = 0, psi(L) = 0, D(psi)(0) = 0, w(0) = 0, w(L) = 0, D(w)(0) = 0, D(w)(L) = 0, u(0) = 0, u(L) = 0;
with(DEtools):
res:=casesplit(sys,[u,w,psi]):
nops([res]);
for i from 2 to 5 do res[i] end do; #Special cases
#############################################################
## Testing the systems at x=0 and with bcs:
for i from 1 to 5 do
  temp:=op(1,res[i]);
  try
    R[i]:=eval[recurse](convert(temp,D),{x=0,bcs});
  catch "numeric exception":
    printf(StringTools:-FormatMessage(cat("i = %d\n",lastexception[2..-1],"\n")),i)
  end try;
end do:
R[2]; # 0 = 27224/3872072675, so no good
R[3]; # 0 = 27224/3872072675, so no good
R[4]; # 0 = 27224/3872072675, so no good
R[5]; # No contradiction
## Looking at the first two equations in res[1]:
temp:=op(1,res[1]):
R[1]:=eval[recurse](convert(temp[1..2],D),{x=0,bcs}); # 0 = 27224/3872072675, so no good
## So far the case res[5] is a possibility (and the only one)
res[5];
## Now look at that case 5 at x=L (and using bcs)
temp:=op(1,res[5]):
eval[recurse](convert(temp,D),{x=L,bcs}); # No contradiction
temp;
## We see that psi is constant (C), so odes for u and w are
sys5:=eval(temp[1..-2],psi(x)=C);
## We solve the w-equation:
bcsw:=op(select(has,{bcs},w));
## Notice that there are 4 boundary conditions
## One of those, w(L)=0, will be used to determine C.
solw:=dsolve({sys5[2],w(0)=0,D(w)(0)=0,D(w)(L)=0}); #Solving exactly
#Now determine C:
eval[recurse](solw,{x=L,w(L)=0});
C_result:=solve(%,{C}); # psi(x) = C
## This is already contradicting psi(L) = 0 which is required in bcs.

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