The PDEtools:-dchange() is made expressly for this purpose.  The attached worksheet shows you how to transform your PDE.  You may need to adjust it a little in order to complete your task.

change-of-variables.mw

It is very likely that internally Maple converts all mathematical expressions to a postfix or prefix form for processing.  I am not aware of a user-accessible hook to that internal representation.

But you will find a very nice external utility for doing that in Notation polonaise inversée.  The part that you would be interested in is the proc toNPI().  For instance, toNPI("a+b*c") returns
    [[identif, "a"], [identif, "b"], [identif, "c"], [binaryoperator, "*"], [binaryoperator, "+"]]
Taking the seond element of each of the five sublists we get
    "a"  "b"  "c"  "*"  "+"
which is close to what you want.

In infix-to-rpn.mw I have extracted the Maple code (but not the comments) from that web page.

I must note that this code represents a very generic "textbook example" of a parser/analyzer in the sense that you would write pretty much the same code in any imperative language such as C or Pascal.  In particular, it does not take advantage of any special features that Maple may offer. I suspect that some parts may be shrunk significantly by calling Maple-specific procedures.

Your assumption that H(x,y)=f(x)*g(y) implies that ln(H(x,y)) = ln(f(x)) + ln(g(y)).  But separating the x and y in the latter form is easy.  This leads to the proc below whose API is:
    ans := doit(H(x,y), [x,y]);
where ans is the list [f(x),g(y)] in case of success or NULL in case of failure.
doit := proc(expr, vars)
  local z1, z2;
  simplify(ln(expr), symbolic);
  expand(%);
  z1, z2 := selectremove(has, %, vars[1]);
  if has(z1, vars[2]) or has(z2, vars[1]) then
    return NULL;
  else
    return simplify([exp(z1), exp(z2)]);
  end if;
end proc:

Examples: (stolen from Kitonum's answer)

doit((3*y + y^2)*3*x/(x + sin(x)), [x,y]);
doit(2^(x^2-y+2*x), [x,y]);
doit((3*y + x^2)*3*x/(x + sin(x)), [x,y]);   # not separable!
doit(s*t*exp(t+s), [t,s]);
doit(sqrt(s*t), [t,s]);
doit(exp(x^2-y), [x,y]);

In particular, sqrt(x*y) does not pose a difficulty with this method.

What is a "circular rod"?  I assume that you mean a rod with a circular cross-section.  That, however, is in conflict with your "1D" requirement.  If it is 1D, how does it matter that the cross section is circular, square, or whatever?

Anyway, if my interpretation of your question is correct, then you are looking at heat conduction in a line segment, let's say 0 < x < L.  The heat equation for the temperature u(x,t) is u_t = u_xx, where the x and t subscripts indicate derivatives with respect to the space x and the time t.  The initial condition is u(x,0)=u₀(x).  You need boundary conditions at x=0, and x=L.  Those depend on the details of what you want to solve.

Your statement "with a heat source at its base" is ambiguous.  Is the temperature of that heat source known?  If so, then the boundary condition at x=0 will be u(0,t)=A(t), where A(t) is the known temperature at x=0.  Alternatively, you may not know the temperature at x=0, but you may know the heat flux, that is, how much heat per unit time is injected into the rod. In that case the boundary condition would be u_x(0,t)=-B(t)., or in Maple notation, D[1](u)(0,t)=-B(t).  Similar considerations apply to the boundary at x=L, but you haven't said anything about that.

This may do what you want:

subs(tmp=1-theta, subs(theta=1-tmp, e1_1));

restart;
with(plottools): with(plots):
Obj := display(dodecahedron([0,0,0], 1));

# PQ is the axis of rotation
P := [1,-1,-1]: Q := [0,2,1/2]:
Axis := pointplot3d([P,Q], connect, color=red, thickness=5):

nframes := 30:
animate(rotate, [Obj, t, [P,Q]], t=0..2*Pi - 2*Pi/(nframes-1),
    frames=nframes, background=Axis,
    scaling=constrained, paraminfo=false);

restart;
de := diff(p(h),h)=A/(B+C*p(h));
ic := p(h0)=p1;
dsolve({de,ic}, p(h), implicit);

Thus, the problem has been reduced to solving a quadratic in the unknown p(h).  If you have some extra information about the problem's constants, then you may pick the desired solution from it.

Answer to question number 1: End each equation with a semicolon.

Answer to question number 2: Write the equations side-by-side, separate them with commas.

Do this:

sol := rsolve({f(n)=f(n-1)+10*f(n-2), f(1)=1, f(2)=2},f(k));
simplify([seq(sol, k=1..10)]);

The result is
[1, 2, 12, 32, 152, 472, 1992, 6712, 26632, 93752]

Solve the PDE numerically:
restart;
pde := diff(u(t,x),t)=diff(u(t,x),x,x);
ic := u(0,x)=sin(x);
bc := D[2](u)(t,-1)=0, D[2](u)(t,1)=0;
pdsol := pdsolve(pde, {ic, bc}, numeric);

Plot the solution at time t=0.2:
pdsol:-plot(t=0.2);

Animate the solution for a sequence of times:
plots:-animate(pdsol:-plot, [t=time], time=0..2,
   font=[Times,20], axesfont=[Default,12], axes=normal);

You should end every statement in Maple with a semicolon.  That includes compound statements with "if", as in:
if a > 1 then
  some statement;
  another statement;
end if;     # note the terminating semicolon of the "if" statement

In certain cases the semicolon may be optionally omitted but you don't have to.  Always terminate your statements with semicolons and you will be alright.

I can't tell what the inner workings of Maple may be on this, but there seems to be some good excuse to justify its behavior.

Maple has no trouble in solving that implicit solution for y(x):
solve(sol, y(x));

However, the structure of the solution is displayed much more clearly if we specify an initial condition:
sol := dsolve({ode, y(0)=1}, y(x));
plot(rhs~([sol]), x=0..4);

We see that the initial value problem is ill-posed—the initial condition does not determine the solution uniquely!  The y(0)=1 is not significant; any initial condition (other than y(0)=0) leads to non-unique solutions.  Thus, asking for a solution y(x) is not a meaningful request, and that's perhaps the reason why Maple stays on the cautious side.

This worksheet mw.mw shows how.  In it I have corrected a typo (a missing negative sign) in your calculations.

That said, it's not clear to me why you bother with Rayleigh-Ritz which was invented way back when, before the age of computers.  A modern numerical solver can obtain a more accurate solution, more easily.

I don't know what the azimuth and polar angles are.  In the attached worksheet I have used the longitude and colatitude angles which are normally used in spherical coordinates.  Perhaps they mean the same thing.

spherical-cap.mw

The problem with your plot command is that it does not refer to sol which contains the numerical solution of the differential equation. What you want is
   
or equivalently
   

You may add time range and color options, as needed.