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