Here is an example:

This example of solving a system of two PDEs is given in Maple's help page on pdsolve/numeric:

PDE := {diff(u(x,t),t)=-1/20*diff(v(x,t),x,x),
        diff(v(x,t),t)=diff(u(x,t),x,x)}:
IBC := {u(x,0)=sin(Pi*x), u(0,t)=0, u(1,t)=0,
        v(x,0)=1-x, v(0,t)=1, v(1,t)=0}:
pds := pdsolve(PDE, IBC, numeric, spacestep=1/50);
plots[display]([seq(pds:-plot(v,t=i/10),i=0..5)]);

Here is an alternative soltution to what vv has suggested.

restart;
with(plots):  with(plottools):
# define the surface z=f(x,y)
f := (x,y) -> sin(x)*sin(y);
surf := plot3d(f(x,y), x=-Pi..Pi, y=-Pi..Pi);
# define an object in space
hoop := tubeplot(1.5*[1/2 + cos(t), sin(t), 2.5 + cos(t)/3],
    t=0..2*Pi, radius=0.1, color=red, style=surface);
# Define the projection map
T := transform((x,y,z)->[x,y,f(x,y)+0.05]):
display([surf, hoop, T(hoop)], scaling=constrained);

That looks like a bug to me, but if we assume mu1 < mu2, we do get an answer in terms of elliptic integrals:

z := 1/sqrt( Pi^2 * gamma1 * gamma2
        *(1+((x-mu1)/gamma1)^2) * (1+((x-mu2)/gamma2)^2) );

int(z, x=-infinity..infinity) assuming gamma1 > 0, gamma2 > 0, mu1 < mu2;

As I commented earlier, it appears that what you have in mind is:

    y[i+1] = y[i] + e(x[i]) mod N
    x[i+1] = x[i] + y[i+1] mod N

Here is a plot of the sequence of 200 iterates, corresponding to N=20 and the initial point at [1,8]:

Here is the worksheet that produced it: mw.mw

Here is the graph of r2 versus s:

The details are in worksheet.mw.

In the worksheet I plot s versus r2.  What I have shown above is r2 versus s which is probably more meaningful. To obtain this alternative graph, change the worksheet's final command to:

plots:-odeplot(dsol, [s, r2(z)], z=-0.94..10,
    labels=['s', 'r2'], numpoints=1000);

You have a system of two second order differential equations, therefore you should expect to have to supply four boundary conditions.  You should not expect, however, to have a solution for arbitrarily assigned values; boundary value problems are  very fussy in that respect.

I assume that you have a reason for being interested in these equations.  The context in which they arise should help you determine the boundary conditions.  In the absence of a context, I picked some arbitrary values for the parameters and four  boundary conditions, and obtained a solution.  Here it is:

nonlinear_ODE-ver2.mw

If this is not to your liking, then you should be more explicit about your needs.

As to the question of "solving faster", we should have some assurance that we are able to solve the system to begin with, and only later worry about making it faster.

This one looks pretty simple to me:

restart;
y := x*sqrt(a*x^2+b*x+c);
int(y, x): collect(%, ln, factor);

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.