In the question p was assumed constant. However, it may be worth mentioning that the non-constant case can be handled as well, but we must use DEtools[mult]:

with(DEtools):
L:=(Dx-p(x))^2;
diffop2de(L,f(x),[Dx,x]); # does not correspond to composition
#This gives the correct version:
mult(Dx-p(x),Dx-p(x),[Dx,x]);
diffop2de(%,f(x),[Dx,x]);

@GOODLUCK If you are rewriting your system using other names, then be sure to avoid gamma, Pi, I, and D, since they have special meanings in Maple:
gamma is Euler's constant, try evalf(gamma);
As said earlier Pi is the usual, try evalf(Pi);
I stands for sqrt(-1), try I^2;
D is the differentiation operator, D(sin); and D(sin)(Pi/3); #compare to diff(sin(x),x);
Also avoid using a name like alpha together with indexed versions of alpha like alpha[0]. Say you assign to alpha[0], alpha[0]:=7.913. Then you have implicitly created a table with name alpha, try eval(alpha);
Also use * as the multiplication sign, the dot is really meant for matrix multiplication.

@GOODLUCK Notice the tilde in my command
fvars:=freeze~(vars);
That means that the procedure freeze will be applied to each element in the set (or list or other "container") vars.
This syntax for elementwise operation was introduced in Maple 13.
See
?elementwise
If instead you just do freeze(vars); then you freeze the whole set of variables.
Notice the different outputs of
freeze(vars); #returns one "freeze" name
and
freeze~(vars); #returns a set of "freeze" names and also regular names (these are not frozen as there is no need for that).

@GOODLUCK fsolve needs to know the values of the parameters in the equations.
I tried this:

params:=indets(feqns,name) minus convert(fvars,set); #The parameters
paramvals:=`=`~(params,.1); #values set arbitrarily to 0.1 each and everyone
#The command
fsolve(eval(feqns,paramvals),convert(subs(Pi=pi,fvars),set));
#didn't finish before my patience (short) ran out.
#Correction: It is worse:
params:=indets(feqns,{name,name^name}) minus convert(fvars,set);
Also you have to be careful about using alpha and alpha[0] (and alpha[whatever] ) in the same calculations.
Evaluating alpha[0] at alpha=.1 :
eval(alpha[0],alpha=.1);
#returns .1[0] which doesn't evaluate to a float.
My suggestion is that you use other names not involving `^` or `[]`.

If possible it is a good idea to treat a considerably smaller system first, the idea being: if a smaller system cannot be solved there is not much chance of success with a larger one.

@GOODLUCK As Carl is pointing out, it is strange to call mu^H a variable. Assuming that by mu^H you don't mean mu raised to the power H, but just mu with an upper index H, you could freeze mu^H so that Maple doesn't treat it as mu raised to power H:

vars:=[c^B,n^B,h^B,mu^B,mu^H,c^s,n^s,h^s,d^s,c^L,m^F,d^L,mu^c,mu^I,
nu^c,nu^I,c^D,d^D,m^R,m^N,K,q,mu^F,k^L,E^L,k^D,E^D,y,M,w^B,P^H,r^B,r^H,
r^D,Pi,r^F,w^F,w^I,w^R,w^N]; #list
fvars:=freeze~(vars); #freezing
feqns:=subs(`=`~(vars,fvars),{eqn||(1..40)}): #Replacing your "variables" by frozen ones
#You used a "variable" called Pi. Pi in Maple is the usual Pi, so use pi:
solve(feqns,subs(Pi=pi,fvars));
No solution.

From the image it looks like there is possibly missing a multiplication sign, so that what is really attempted is
EC2:= Z=1+beta-(q*beta*(Z-beta)/(Z+epsilon*beta)(Z+sigma*beta));
solve(eval(EC2, P), Z);

In Maple 15 this results in the warning mentioned, whereas in Maple 16 and 17 it results in errors (different though).

From the image it looks like there is possibly missing a multiplication sign, so that what is really attempted is
EC2:= Z=1+beta-(q*beta*(Z-beta)/(Z+epsilon*beta)(Z+sigma*beta));
solve(eval(EC2, P), Z);

In Maple 15 this results in the warning mentioned, whereas in Maple 16 and 17 it results in errors (different though).

Why don't you just write the ode in Maple syntax and use dsolve?

@Axel Vogt Thank you!
You forgot "return" before 'procname'.
I still wonder what goes on in plot and plot3d.

@Axel Vogt Thank you!
You forgot "return" before 'procname'.
I still wonder what goes on in plot and plot3d.

@Axel Vogt There is a difference in the way plot and plot3d handles modp:
restart;
modp(5.4321,2); #An error as expected
plot(x->modp(x,2),-3..3); #Works, surprisingly
plot('modp(x,2)',x=-3..3); #Works, surprisingly
plot3d((x,y)->modp(x,2),-3..3,-1..1); #Doesn't work
plot3d('modp(x,2)',x=-3..3,y=-1..1); #Works

@Axel Vogt There is a difference in the way plot and plot3d handles modp:
restart;
modp(5.4321,2); #An error as expected
plot(x->modp(x,2),-3..3); #Works, surprisingly
plot('modp(x,2)',x=-3..3); #Works, surprisingly
plot3d((x,y)->modp(x,2),-3..3,-1..1); #Doesn't work
plot3d('modp(x,2)',x=-3..3,y=-1..1); #Works

Your solution works, but I wonder why since
F(5.78);
results in the error

Error, (in F) invalid argument for modp or mods

as does also
modp(5.78,2);

Your solution works, but I wonder why since
F(5.78);
results in the error

Error, (in F) invalid argument for modp or mods

as does also
modp(5.78,2);

I fail to understand why you expect answers of type numeric in what you are doing. In fact I don't understand what you are trying to do.
Your equations eq1 and eq2 could be obtained in a somewhat more readable form like this:

restart;
egf := exp(2*x*z-z^2);
test1r := subs(x=dummy,z=x,egf);
test1r2 := solve(test1r=f,x);
b := Matrix([[a1,a2],[a3,a4]]);
invb := b^(-1);
testing1:=<f|f^2>.invb.Matrix([[x,1],[1,x^2]]).b; #Row vector
eq1 := test1r2[1] = testing1[1];
eq2 := test1r2[1] = testing1[2];