The series
S:=4/Pi*Sum((-1)^(n+1)/(2*n-1)*cos(2*Pi*(2*n-1)*x),n=1..infinity);
is a fourier series of a simple function g.
This function is even and has period 1. In fact on the interval 0..1 it is given by

g:=piecewise(x<1/4,1,x<3/4,-1,1);

These facts are readily verified the usual way.
Thus in h3 you can replace the series S above by g:

H3:=1+lambda*m*x/a0+phi*4/Pi*Sum((-1)^(n+1)/(2*n-1)*cos(2*Pi*(2*n-1)*x),n=1..infinity);
h3:=subs(Pi/4*S=Pi/4*g,H3); #S is not an operand of H3, but Pi/4*S is!
Having done that and after having assigned values to the parameters you can just integrate:
DP3val:=value(DP3);
which is obviously easy to plot.

INT:=Int(1/sqrt(a*x^3+1),x=0..X);

On a range in which a*x^3+1 is positive INT is an increasing function of X whether or not X<0 or X>0.
Even when you start with a concrete value for a the answer is wrong for X>0:

value(eval(INT,a=0.1)) assuming X>0;
plot(%,X=0..1); #Decreasing
#This is OK though:
value(eval(INT,a=0.1)) assuming X<0;
plot(%,X=-1..0);

Try replacing the last few lines by

params:={a=1,b=1,c=1};
 f:=(x,y)->x*(x-1)*(y-1);
 g:=(x,y)->0;
uc:=(x,y,t)->eval(u(x,y,t),params);
U:=(x,y,t)->Sum(Sum(uc(x,y,t),n=1..20),m=1..20);
plot3d(U(x,y,2),x=0..1,y=0..1,axes=boxed);

#Specifically: don't assign to the variables you assumed on:
restart;
assume(n,integer);
u:=n;
n:=78;
u;

N:=10 ; # Global Var
F:=(x,y)->signum(abs(x-N/2)+abs(y-N/2)-N/4); #signum
Average := proc(F, f0::Matrix) local f, i, j;
  f := copy(f0); # copy
  for i to N do
    for j to N do
      if F(i, j) < 0 then
        f[i, j] := (f0[i - 1, j] + f0[i + 1, j] + f0[i, j + 1] + f0[i, j - 1])/4 ;
      end if
    end do
   end do;
   f
end proc;
f0:=Matrix(N,F);
Average(F,f0);
f0;

Have you had a look at convert/confrac?

Here are two different methods:

restart;
eq:=58.50208000*omega^2-5.263157894*10^11-3.150112000*10^9*Pi^2+5.703952800*a*omega^4-5.131578947*10^10*a*omega^2-4.033564463*10^8*a*omega^2*Pi^2-8.638197750*10^7*a^3*Pi^4*omega^2+7.771381575*10^17*a^3*Pi^4+4.651337250*10^15*a^3*Pi^6+5.181105261*10^15*a*Pi^4 = 0;
res:=solve(eq,a); #solving for a: 3 solutions (cubic in a)
plots:-complexplot(res[1],omega=0..5,thickness=3); #real values as it appears
plots:-complexplot(res[2],omega=0..5,thickness=3);
plots:-complexplot(res[3],omega=0..5,thickness=3);
plot(res[1],omega=0..5000,0..1e-5);
#Using implicitplot instead
plots:-implicitplot(eq,omega=0..5000,a=0..1e-5);

The error that you get tells the story. There are too few boundary conditions. Your equation is third order in X, so you need 3 besides the one initial condition.
You can try this, where I rather arbitrarily added a boundary condition:

KDV2:= diff( u(X,T), T)+ 6*u(X,T)*diff(u(X,T),X)+ diff(u(X,T),X$3);
inf:=4:
SOL:=pdsolve(KDV2,{u(-inf, T) = 0, u(inf, T) = 0, D[1](u)(inf,T)=0, u(X, 0) =1},numeric,time=T,range=-inf..inf,spacestep=inf/50);
SOL:-animate(T=0..1);


Notice that infinity is set to 4 instead of 10. I had no immediate luck with 10.

If your constant Lambda is positive then there is exactly one real root.

restart;
u:=F+(1/4)*(4*(-Lambda)^(3/2)*h*m^2-2*sqrt(-Lambda)*exp(Lambda*h^2*m^2)*h*m+sqrt(Pi)*erf(h*m*sqrt(-Lambda)))/(m^2*(-Lambda)^(3/2));
us:=simplify(u) assuming Lambda>0;
limit(us,m=0); #The limit exists and is h+F
##We claim that the derivative is positive for all m if h>0 (and negative if h<0):
d:=simplify(diff(us,m));
convert(d,dawson);
factor(%);
#####
Let us assume that we know that x^3-x+dawson(x)  (which is odd) is positive for all x>0.
To show that requires little work, but just take it for granted now. (But see below).
Then we see that if h>0 then the derivative d is positive for all m since it is even in m. It tends to infinity as m-> infinity. Thus u is increasing in m and has range -infinity..infinity.
If h<0 instead similarly we find u to be decreasing with the same range.
Thus u has exactly one root.
To find the root use fsolve.
Illustration of the fact used about dawson:
plot(x^3-x+dawson(x),x=-1..1);
Proof: For x>=0 we have dawson(x)=exp(-x^2)*Int(exp(t^2),t=0..x) >= exp(-x^2)*Int(1,t=0..x) =
   exp(-x^2)*x >= (1-x^2)*x = x-x^3.
Thus x^3-x+dawson(x) >=0 for all x>=0. In fact equality only if x=0.

You need initial conditions as Carl says, and you need to know T0 since you must solve numerically.
Here is an example:

res:=dsolve(eval({EQ1,EQ2,EQ3,EQ4,EQ5},T0=1) union {q1(0)=0.1,q2(0)=0.2,q3(0)=0.3,q4(0)=0.1,q5(0)=0.2,D(q1)(0)=0,D(q2)(0)=0,D(q3)(0)=0},numeric,maxfun=10^6);
plots:-odeplot(res,[seq([t,cat(q,i)(t)],i=1..5)],0..5e-4);

Maybe I misunderstand your problem, but does this not do what you want:

restart;
Array1 := Array([Array([1, 2]), Array([3, 4])]);
Array2:=copy(Array1);
Array1[1]:=Array([11,45]);
Array1;
Array2;

I think I see your reason for using vector in two cases as opposed to Vector. However, using vector should be avoided as it is the old structure used in the deprecated linalg package.
You are using the vector elements as names only, so there is no need for vectors or Vectors. Just use indexed names.
If you do insist, you can use Vector with symbol=some name as shown in the two lines that have been commented out.

restart;
m0 := t-> 1-exp(-t*.5);
m := t->(1/(1+exp(-t+5))-0.67e-2)*1.0067;
n[max] := 10; delt := .2; n := n[max]/delt;
T := Vector(50);
#b := Vector(50,symbol=bb); #No need for this
for i from 2 to n do T[i] := T[i-1]+delt end do:
fun := t->add(b[i]*m0(t-T[i]), i = 1 .. n);
#fu := Vector(50,symbol=f); #No need for this
for x to 49 do fu[x] := fun(x*delt) = m(x*delt) end do:
s := solve({seq(fu[i],i=1..5)}, {seq(b[i],i=1..5)});

The basic reason is that in general the simplification would be wrong. Specifically it is wrong whenever L<=0.

Try these:

ineq:=L*xi<=L/2;
solve(%,xi) assuming L>0;
map(x->x/L,ineq); #You take responsibility for the soundness of the operation
isolate(ineq,xi) assuming L>0;

Does this idea work:

restart;
with(IntegrationTools);
Z:=a*Int(f(x),x=0..1);
Combine(Z); #a stays outside
Z0:=Int(b,x=0..1); #b unassigned
eval(Combine(Z0+Z),b=0);
GetIntegrand(%);
Z:=-a*Int(-f(x),x=0..1);
eval(Combine(Z0+Z),b=0);
GetIntegrand(%);

This seems to me a bug, as also hinted at in my comment to the question.

restart;
A := 54.15836673*(diff(X2s(x), x, x)) = -365.4395362*(diff(X1s(x), x))+208.2315661*X2s(x);
B := 641.1196154*(diff(X1s(x), x, x)) = 365.4395362*(diff(X2s(x), x))-2.575699975*X1s(x)-7.882173342;
indets(dsolve({A,B}),name) minus {x}; #Output {_C1, _C2, _C3}
sysE:=convert([A,B],rational,exact);
indets(dsolve(sysE),name) minus {x}; #Output {_C1, _C2, _C3}
indets(dsolve({A,B},convert_to_exact=false),name) minus {x}; #Output {_C1, _C2, _C3,_C4}
A1 := a*(diff(X2s(x), x, x)) = -b*(diff(X1s(x), x))+c*X2s(x);
B1 := d*(diff(X1s(x), x, x)) = b*(diff(X2s(x), x))-e*X1s(x)-f;
indets(dsolve({A1,B1}),name) minus indets({A1,B1},name); #Output {_C1, _C2, _C3,_C4}

I would have thought that the correct output {_C1, _C2, _C3,_C4} would at least come in the exact cases, but possibly not in the convert_to_exact=false case.

It should be noted that
dsolve({A,B},{X1s(x),X2s(x)},method=laplace):
{X1s(0),X2s(0),D(X1s)(0),D(X2s)(0)} subset indets(%,function);
returns true. Thus also here there are 4 constants.

I shall submit an SCR.

The singularity warning means here that x(t) and/or y(t) reaches infinity in finite time.
With inits being defined as your list of initial conditions you can do

DEplot(DE13, [x(t), y(t)], t = -5 .. 2, inits, x=-5..5,y=-2..2,
stepsize = 0.1e-1, title = "phaseplane 3 prime plot", linecolor = black, thickness = 1);