I am a noob in (P)DE. Anyway: divide by F5 (not zero?) to get rid of it and just consider
F3/F5 as part of phi4 - which is legitimate here.

Now you have 0 = pde(r,t) + arbitraryFunction(r,t). Or otherwise said: pde = arbitraryFunction(r,t).

How should that give an explicit solution?

In your 2nd case this is not arbitrary, it is a pde for phi4

One always can give a 'trivial' solution, the polynomial of degree n-1 given by
that n pairs (j, List[j]), j = 1 .. n. It has integer coefficients in any case.

g(x) =  -487+488*x-485*binomial(x,2)+487*binomial(x,3)-490*binomial(x,4)+
             488*binomial(x,5)-457*binomial(x,6)+337*binomial(x,7)

But that's not what you have in mind, g(8) = 37, g(9) = 845, g(10) = 5554

Have you tried to save a sheet as HTML? Where I sometimes got better result by importing graphics through Open Office Writer instead of MS Word. Or mark the graphics with the mouse and right click to get a context menu to save it in a convenient format

I never work with Units, but:

restart;
occultRatio := proc(beta, R, alt)
 local rslt, cb, rr; 
  rslt:=0;
  cb:=cos(beta);
  rr:=sqrt(1-(R/(R+alt))^2);
  if (cb > rr) then
   rslt := abs(-Pi/2 + arcsin( rr/cb ))/Pi;
  end if;
  return rslt;
end proc;

JupiterRadius := 69173;
OrbitRadius := 2* JupiterRadius;

plot('occultRatio'(b, JupiterRadius, JupiterRadius), b = 0 .. Pi/2);
plot('occultRatio'(b, 2, 2), b = 0 .. Pi/2);

I do not want to copy and clean up your code, but here is a suggestion:

Do not use floats like rho := 0.122e4, but append rho:=convert(rho, rational)
and try agian

K[x] / (x^2) should be a canonical way for 'dual' numbers (as Wiki) mentions

Not sure for quaternions (never needed it being off Physics)

What are m and phi in your 2nd case? If it has no numerical value or a reasonable
property, then evalf can not know about the imaginary part.

Remark1: there is no need to use 'VectorCalculus' and the complicated input here,
you can just type in the formulae. This is just evalf(t1a) and then picking those with
positibe imaginary part.

Remark2: in the 1st case beta,Q,P are not used and not needed, dito in the 2nd
case - but you may wish to provide m and phi

Remark3: setting phi = psi/m will make it more simple

Have you looked up the help, ?Ei

So I guess your question is: plot the zeroes of f(r,s, t) for r,s in a given range.

This will be a surface (depending on r,s).

doodles.zip

This is more or less something exponential, so log may be the way
for a reasonable scaling, up to signs

plot([ln(temp6c),ln(-temp6c)],Sigma=-40..40);

log(zero) = -infinity, thus it is peaked where you may guess some zero

Here is a sketch, how I try to understand pagan's solution:

If we agree to reduce to the (open) unit hypercube and that max( ()^2 ) = max(..)^2
then I think one can argue as follows:

The set, where any of the coordinates may coincide, is low dimensional (analytic of
dim <= n-1) and thus does not contribute to the integral.

Denote the complement by U.

Let be U0 := points with 0 < x1 < x2 ... < xn < 1. The symmetric group Sn operates by
permutation on the coordinates. Any of those maps U0 to U. Any point in U has just one
ordering of its coordinates, thus is an image of a point in U0 by a unique permutation.

Thus U is the disjoined union of images of U0 under the symmetric group.

And the integral is the sum over those images.

Now we have to convince ourselfs, that this works reasonable on the integrals.

The permutations do not change volume. And on the function max(...) they fit, since
max is invariant under permutations. Therefore the integrals in the sum have all the
same value, that over U0.

There are as many summands as elements of Sn, hence n!-times that over U0.

The integral over U0 is just as already written down by pagan for n=5, 0< v<w<z<y<x < 1
as it is seen by 'resolving the range' on paper.

That is a bit of handwaving, but probably can be made sound.

see the attached file hirnyk_max_squar.mws

roughly -0.19138848e-3 + 0.60440364e-4*I

First: if you allow the lower bound to be variable, then you may get a solution
for nay of those bounds, i.e. a curve. Not sure, whether that is intended.

Here is a way to do it for b=0 and I prefer to use rational constants, not floats.
And please note the UpperCase for the integral nad evalf on it (you may wish to
read that topics in the help, it is worth to do):

restart;
Digits:=15;

g:= proc(u) 
  evalf( Int((1-x)*(1-(1/5*u*(1-u)/x/(1-x)-1/5)^(1/4)),x = 0 .. u, method = _d01ajc) ) - 1/20  
end proc;

plot(%, u = 0 .. 1/2, numpoints=10);

You want a zero of g and there is one close to u = 0.2

For that I use my own solver (Brent) and get u = 0.197213988275156 as solution.

May be you want to plot the coordinates of the points in the plane?

L := ...
M:=map(convert,%, list);

plot(M, style=point, symbolsize=20);
plots[pointplot](M, symbolsize=20, color = blue);

In case you want to work with Vectors you may want to spent some
time reading the help, how they can be used,  <0,0> is more easy
than your notation.