Using method= _MonteCarlo and given a large enough epsilon:

int(add(cat(x,i),i=1..10)^2,seq(cat(x,i)=0..1,i=1..10)],numeric,method=_MonteCarlo,epsilon=0.5e-4);

I got 25.83247877.

First an observation: In Maple 8 the 3 results of the evalf command in the following have no imaginary parts at all:
 

M:=MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi);
seq(evalf(eval(M,x=xx)),xx=[0.4,0.5,0.6]);

In Maple 12, 15, 16, ..., 2018 there are very small imaginary parts in all 3.
This explains the difference in the results of this workaround for the integral:

u:=(-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2);
ode:=diff(y(x),x)=u;
res:=dsolve({ode,y(4/10)=0},numeric); #OK in Maple 8, but not in 12, 15...2018
##A workaround for the workaround in 12, 15...2018
res:=dsolve({ode,y(4/10)=1e-22*I},numeric); 
res(6/10);

The result of res(6/10) in Maple 2018.1 is -0.555168575558567e-2+1.00000000000001*10^(-22)*I
with the default settings of abserr and relerr.
 

I tried the following and lost the kernel connection.
 

restart;
u:=(-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2);
evalf(u);
int(%, x = 4/10 .. 6/10, numeric);

I shall report that as an SCR.

I once had that problem and then made a small procedure I called Chop lacking a better name.
Unlike evalf[n] which will turn exact expressions like Pi or sqrt(2) into floats with n digits Chop only works on existing floats.
 

Chop:=proc (x::anything, n::posint) local op1, op2, k; 
   if not hastype(x, float) then return x end if; 
   if type(x, float) then 
      op1, op2 := op(x); 
      k := length(op1); 
      evalf(round(op1*10^(-k+n))*10^(op2-n+k), n) 
   elif type(x, complex(float)) then 
      procname(Re(x), n)+I*procname(Im(x), n) 
   else 
      evalindets(x,float,procname,n) 
   end if 
end proc;

Examples:
 

Chop(123.456789,4);
Chop(1.23456789,4);
evalf[10](f^Pi);
evalf[4](%); #Actually doesn't do anything!
Chop(%,4);

I made a .txt file and saved it at F:/test.txt.
Then I tried the example in the help page for readline, but replaced "infile.text" with the full path to my file.
Then no problem:
 

restart;
file:="F:/test.txt";
first := readline(file);
line := readline(file);
while line <> 0 do
   last := line;
   line := readline(file)
end do;

Try this:
 

f:=diff(ChebyshevT(n, r), r);
g:=simplify(eval(f, r = 0));
simplify(g) assuming n::odd;
simplify(eval(%,n=2*k+1)) assuming k::integer;
simplify(g) assuming n::even;

@9009134 The numerical solver pdsolve/numeric can solve a system of pdes that in vector form can be written as
diff(u(x,t), t) = f(t, u(x,t), diff(u(x,t),x), .., diff(u(x,t),x$n))
where u and f are vector valued.
If the given system is not first order in t it is converted to one (see the help page pdsolve/numeric).
Initial conditions of the form u(x,t0)=g(x) and boundary conditions must be given on two lines x=a and x=b.
The problems handled are thus time-based, i.e. the algorithm starts with t=t0 and u(x,t) = g(x) and then will work its way from there to any given t and along the way making sure that the boundary conditions are satisfied.

Your problem is not formulated that way. You have conditions given on the rectangle given by the lines x=0, x=L, z=0, z=h.
Your system could be rewritten (by pdsolve itself) as a system of the form given above if x is considered time, but then you cannot impose conditions on both of the lines x=0 and x=L. You must choose one of them, say x=0.  Furthermore since your system is second order in x (time) in both p1,p3, and Phi you must give initial conditions on the time derivatives too for x = 0.
A one dimensional analogue of the difference between your pde-bvp problem and the time based problem described here is the difference between a bvp problem for an ode and an initial value problem for the same ode.

### To make my comment clearer let me try doing what I said, i.e. consider x time and use z as the space variable. Notice my change of pd3. You probably meant one of the two occurrences of diff(Phi(x,z),x,x) to be diff(Phi(x,z),z,z). (If you meant what you wrote you must remove one of the boundary conditions for Phi).

restart;
pd1 := -alpha*beta^2*diff(p1(x, z), x, x)-alpha*beta1^2*diff(p3(x, z), x, z)-alpha*beta2^2*diff(p1(x, z), z, z)-alpha*beta3*diff(p3(x, z), x, z)+alpha*p1(x, z)+diff(Phi(x, z), z) = f1(x)+z*f2(x);
pd2:=-alpha*beta1^2*diff(p1(x, z), x, z)-alpha*beta^2*diff(p3(x, z), z, z)-alpha*beta2^2*diff(p3(x, z), x, x)-alpha*beta3^2*diff(p1(x, z), x, z)+alpha*p3(x, z)+diff(Phi(x, z), z) = g(x);
##
## I have corrected your equation pd3:
## I changed one of the two occurrences of diff(Phi(x,z),x,x) to diff(Phi(x,z),z,z).
##
pd3 := -xi*(diff(Phi(x, z), x, x)+diff(Phi(x, z), z, z))+diff(p1(x, z), x)+diff(p3(x, z), z) = 0; 
alpha := 1; beta1 := 2; beta2 := 3; beta3 := 4; h := 1; L := 5; xi := 6; beta := 4; 
f1 := x->1/cos(x)^2; 
f2 := x-> 1/(sin(x)^2+x^4); 
g := x-> 1/cos(x)^4;
bcs:= Phi(x, 0) = 0,Phi(x, h) = 0, p1(x, 0) = 0, p1(x, h) = 0, p3(x, 0) = 0, p3(x, h) = 0; 
ics:= Phi(0, z) = 0, p1(0, z) = 0, p3(0, z) = 0, D[1](Phi)(0,z)=0, D[1](p1)(0,z)=0, D[1](p3)(0,z)=0 ;
res:=pdsolve({pd1,pd2,pd3},{bcs, ics},numeric);
res:-plot3d(Phi,x=0..0.1);

Maybe it is helpful to show the simple conversion of your system to a system that it is of the form
diff(u(x,t), t) = f(t, u(x,t), diff(u(x,t),x), .., diff(u(x,t),x$n)).
 

restart;
pd1 := -alpha*beta^2*diff(p1(x, z), x, x)-alpha*beta1^2*diff(p3(x, z), x, z)-alpha*beta2^2*diff(p1(x, z), z, z)-alpha*beta3*diff(p3(x, z), x, z)+alpha*p1(x, z)+diff(Phi(x, z), z) = f1(x)+z*f2(x);
pd2:=-alpha*beta1^2*diff(p1(x, z), x, z)-alpha*beta^2*diff(p3(x, z), z, z)-alpha*beta2^2*diff(p3(x, z), x, x)-alpha*beta3^2*diff(p1(x, z), x, z)+alpha*p3(x, z)+diff(Phi(x, z), z) = g(x);
pd3 := -xi*(diff(Phi(x, z), x, x)+diff(Phi(x, z), z, z))+diff(p1(x, z), x)+diff(p3(x, z), z) = 0;
### Now do:
S1:={diff(p1(x,z),x)=q1(x,z),diff(p3(x,z),x)=q3(x,z),diff(Phi(x,z),x)=F(x,z)};
S2:=subs(S1,{pd1,pd2,pd3});
S2i:=solve(S2,{diff(q1(x,z),x),diff(q3(x,z),x),diff(F(x,z),x)});
### The system consisting of S1 and S2i together is now on the desired form:
map(print,S1 union S2i);

Some of your equations in eqs are sets. If you replace eqs by EQS as defined below you will get another error:

Error, (in dsolve/numeric/DAE/make_proc) number of unknown functions and equations must match, got 20 functions {C[0], C[1], C[2], C[3], C[4], C[5], Tg[0], Tg[1], Tg[2], Tg[3], Tg[4], Tg[5], Tg[6], Ts[0], Ts[1], Ts[2], Ts[3], Ts[4], Ts[5], Ts[6]}, and 19 equations
 

## Use this after your final definition of eqs:
EQS:=evalindets(eqs,set,op); # this removes the curly brackets
## Now use
sol := dsolve([op(EQS), op(ICs)], numeric, range = 0 .. tmax, abserr = 10^(-3), relerr = 10^(-3), stiff = true, maxfun = 0);

You must now think about the error message quoted.

The method as well as the notation used seem to be inspired by a Maple webinar by Sam Dao available from Maplesoft: Recorded webinars.

Your odes are linear and the boundary value problems can probably all be solved exactly by dsolve in Maple directly. So these are obviously intended as test cases.
Kitonum solves the initial value problem u(0)=0, D(u)(0)=a exactly and use solve to find the exact value of a by requiring u(1)=0.
I shall give a numerical solution by a shooting method. 
Here I shall use Problem 6:
 

restart;
ode:=(1-x^2)*diff(u(x),x,x)-4*x*diff(u(x),x)-(1+x^2)*u(x)=x^2+cos(3*x);
sol:=dsolve({ode,u(-4)=8,u(-2)=-1}); #Exact solution directly
plot(rhs(sol),x=-4..-2);
## Now we find the numerical solution by shooting:
res:=dsolve({ode,u(-4)=8,D(u)(-4)=u1},numeric,parameters=[u1],output=listprocedure,abserr=1e-10,relerr=1e-9);
Up:=subs(res,u(x));
Q:=proc(u1) res(parameters=[u1]); Up(-2)+1 end proc;
plot(Q,0..10);
s:=fsolve(Q,0..20);
res(parameters=[s]); #Set the parameter u1 at s
plots:-odeplot(res,[x,u(x)-rhs(sol)],-4..-2); # The error

Finally it should be said that dsolve/numeric can solve the boundary value problem without using the shooting method. Just use
 

res2:=dsolve({ode,u(-4)=8,u(-2)=-1},numeric);

######################################################
To get tabular results you can do as follows, where I use the exact solution sol and the approximate solution res (and concretely Up) found above.
 

U:=unapply(rhs(sol),x); # The formula for the exact solution
N:=10: # A choice
interface(rtablesize=N+2); 
X:=Vector([seq(-4+i*2/N,i=0..N)],datatype=float); # The x-values
Uexact:=evalf[15](U~(X)); # The exact solution at those x-values
Uapprox:=Up~(X); # The approximate solution  at those x-values
Udiff:=Uapprox-~Uexact; # The error: Approx-Exact
A:=<X|Uapprox|Uexact|Udiff>; # The data 
L:=Array([x,u_approx(x),u_exact(x),u_approx-u_exact]); # Description
resM:=Matrix([[L],[A]]); # Final result

I don't have Maple 7 on this machine, but I do have Maple 8.
I see all 4 lines in all 3 displays.
Since it appears that things are different in Maple 7 you could try using a list rather than a set. A list respects the order.
Even a sequence should work (it does in Maple 8).
It may be a matter of one plot covering the other?
So try permutations of:
 

plots[display]( [lin[1],lin[2],lin[3],lin[4], squ,poinC,tetA], scaling=constrained,axes=normal);

A slightly modified version of your code that works in Maple 2018.1 too:
 

restart:
 # # # # # # # # # # # # # # # # # # # # # # # # # # # #
 # Areas of portions of a square 
 # # # # # # # # # # # # # # # # # # # # # # # # # # # #
 with(plots):
 with(plottools);
 side:=6:A1:=16:A2:=20:A3:=32:
 # Pick a point in the square: say (x1, y1)
 poinC:=point([side/4,side/3]):
 
poinC:=point([side/4,side/3]);
cpoin:=convert(op(1,poinC),list);
# x1:=cpoin[1];
# y1:=cpoin[2];
 lin[1]:=line([side/2,0],cpoin, color=black):
 lin[2]:=line([0,side/2],cpoin, color=black):
 lin[3]:=line([side/2,side],cpoin, color=black):
 lin[4]:=line([side,side/2],cpoin, color=black):
 #textA:=textplot([(cpoin[1]+side/2)/2,1,`A1`],align={ABOVE,RIGHT}):
 #textA:=textplot([(1.5+side/2)/2,1,`A1`]):
 tetA:=textplot([1,1.3,`A1`]):
 
 squ:= polygon([[0,0], [0,side], [side,side],[side,0]], color=white):
 # Plot without tetA
 plots[display]({seq(lin[i],i=1..4),squ,poinC}, scaling=constrained,axes=normal);
 # Plot with tetA
 plots[display]({seq(lin[i],i=1..4),squ,tetA,poinC}, scaling=constrained,axes=normal);
 # Plot with tetA (placed  differently)
 plots[display]({squ,poinC,tetA,lin[1],lin[2],lin[3],lin[4]}, scaling=constrained,axes=normal);

pdsolve/numeric even insists that boundary conditions must be given on a region defined by lines parallel to the axes. Thus you cannot have a boundary condition like sh(H(t), t) = 0, unless H(t) is a constant, i.e. even if H(t) is actually known (e.g. H(t)=t), it wouldn't be allowed.

I cannot tell you what goes wrong. The logarithmic terms could possibly become imaginary at some point?
Anyway your pde1 can be considered an ode with the parameters f1(x), f2(x), f3(x), f4(x) as is done below.
The ugly and weird notation `σh3` and `σv3` is the 1D rendering of your 2D code. I hate 2D input with a passion. You see why.

pde1;
ode:=subs(Lambda(x,t)=Lambda(t),f1(x)=F1,f2(x)=F2,f3(x)=F3,f4(x)=F4,pde1);
res:=dsolve({ode,Lambda(1.8936825887327868318*10^15) = F4},numeric,parameters=[F1,F2,F3,F4]);
## Testing. First setting parameters. We use x = 1000 here:
res(parameters=[F1=`σh3`(1000),F2=`σ`(1000),F3=Jir3(1000),F4=Lambda3(1000)]);
res(Tf); #result
## A procedure to do that:
LL:=proc(xx,tt) if not type([xx,tt],list(realcons)) then return 'procname(_passed)' end if;
  res(parameters=[F1=`σh3`(xx),F2=`σv3`(xx),F3=Jir3(xx),F4=Lambda3(xx)]);
  subs(res(tt),Lambda(t))
end proc;
## Testing the procedure:
LL(1000,Tf);
## Plotting
plot(LL(1000,t),t=Ts..Tf); # x =1000
plot(LL(x,(Ts+Tf)/2),x=Xvp..Xp); # t= (Ts+Tf)/2
## The animation takes a while and is not very illuminating as it is:
plots:-animate(plot,[LL(x,t),x=Xvp..Xp],t=Ts..Tf);

Surely this couldn't be the most efficient way of producing a solution.

Yes, I give you a very simple example where the initial and boundary conditions are defined by numerical procedures from dsolve/numeric:
 

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x,x);
ibcs:={u(x,0)=f(x),u(0,t)=g1(t),u(1,t)=g2(t)};
#Test with simple conditions
test_res:=pdsolve(pde,eval(ibcs,{f(x)=x,g1(t)=0,g2(t)=t+1}),numeric);
test_res:-animate(t=3);
# Now an odesys producing in fact the same functions f, g1, and g2:
odesys:={diff(f(x),x)=1,diff(g1(x),x)=0,diff(g2(x),x)=1,f(0)=0,g1(0)=0,g2(0)=1};
oderes:=dsolve(odesys,numeric,output=listprocedure);
F,G1,G2:=op(subs(oderes,[f(x),g1(x),g2(x)]));
res:=pdsolve(pde,eval(ibcs,{f=F,g1=G1,g2=G2}),numeric);
res:-animate(t=3);
test_res:-value()(.1234,3);
res:-value()(.1234,3);

Yes, this is a bug. The "solution" found by Maple 2018.1 (and also by Maple 2017.3) is
sol := u(t, x, y) = exp(-t*(x^2+2*x*y+y^2+1))
but it is wrong. It satisfies the initial condition, but not the pde.

In earlier versions (Maple 2016 and earlier) no solution was found.

I shall submit a bug report (aka an SCR).

This is clearly a bug in pdsolve/numeric.
A workaround for your case is to use the following pde (PDE2), which actually on the given range is identical to yours:

PDE2:=PDE+(0=diff(L1(x,t),x)*Heaviside(-x));

If you try

simplify(PDE2) assuming x>0;

you get PDE, but don't use that simplified version, keep PDE2 and use IBC := {L1(Xp, t) = 1, L1(x, Ts) = 1}.
Here is an extremely simple example of this same bug:
 

restart;
pde:=diff(u(x,t),t)=u(x,t);
ibcs:={u(x,0)=x,u(0,t)=0};
pdsolve(pde,ibcs,numeric,time=t,range=0..1);
ibcs2:={u(x,0)=x};
pdsolve(pde,ibcs2,numeric,time=t,range=0..1);
sol:=pdsolve({pde,u(x,0)=x}); #Exact version
##### The trick:
pde2:=pde+(0=Heaviside(-x-1)*diff(u(x,t),x));
simplify(pde2) assuming x>=0;
ibcs:={u(x,0)=x,u(0,t)=0};
res:=pdsolve(pde2,ibcs,numeric,time=t,range=0..1);
res:-plot(t=5);
res:-animate(t=0..5);
## The same for the exact solution
plots:-animate(plot,[rhs(sol),x=0..1],t=0..5);

I think I may have sent an SCR (bug report) a long time ago about this, but I shall do it again just in case.