Is there a way to solve this using pdepe?
No.
Don't you need a total of 4 boundary conditions instead of only one ?
Are you sure about a mixed-term d^2(phi_s)/dt dz in the boundary condition ?
Does xi(u_dot+2alpha_dotz) and G(u-alphah0) mean xi*(...) and G*(...) or xi@(...) and G@(...) ?

