I don't know how to rebind keys, but Ctrl-End takes you to the end of the worksheet on linux.

As I said my previous post, obtaining the optimal path of the boat is a calculus of variations problem.  It turns out that although the problem is doable in principle, the calculations would be quite difficult without the help of a symbolic computational software such as Maple.  This is a good illustration of the non-trivial use of a modern CAS.

The problem's statement and solution are given in detail in the attached worksheet.  A brief version is this.  We have a river represented by the strip 0 < x < 1,   -∞ < y < ∞ in the Cartesian x-y plane.  The water velocity is <0, V(x)>.  A boat travels from (0,0) to (1,0) with a constant speed S relative to the water.  What path to take to minimize travel time?

Here is the shape of the optimal path corresponding to V(x) = 3*x*(1-x) and S=1.  The travel time turns out to be T=1.15.

Here is the worksheet: boat2.mw

Added later: I found and corrected a few typos in the previous worksheet's comments (not the math). Here is the corrected version: boat2-ver2.mw

Earl, now I am able to read your worksheet.  I haven't examined it carefully because while I was waiting for it, I made a worksheet of my own which appears to be quite similar to yours.  Here it is:

boat.mw

In this worksheet the river is the vertical strip 0 < x < 1 in the x-y Cartesian plane and the water flows in the positive y direction with a velocity which is a linear function of the distance from its mid point to either bank, although any other velocity profile may be substituted for it.  The worksheet has a more detailed description of the problem and its assumptions.

The boat goes from (1,0) to (0,0) while pointing toward the destination (0,0) at all times. The following animation, extracted from that worksheet, shows the boat's path:

This answers only one of your two questions.  The other is to find the fastest path between the two destinations. That appears to be quite closely related to the brachistochrone problem from the calculus of variations.  If I find some time, I will work it out and post, unless someone else with more free time on their hands beats me to it.

The code you have posted is not quite readable.  You wouldn't want to read someone else's code that looks like that, would you?

Try to help people who you want to help you.  Format your code as best as you can when you post it.

OK, I did go through the trouble of formatting your code, and simplified it.  Here is what it does.

restart;
pde := diff(u(x,y,z,t),t) + 6*u(x,y,z,t)*diff(u(x,y,z,t),x)
       + diff(u(x,y,z,t),x,y,z) = 0;

Apply the following change of variables:

tr := {t = T/a^beta, x = X/a^alpha, y = Y/a^mu, z = Z/a^nu,
       u(x, y, z, t) = U(X, Y, Z, T)/a^Zeta};
PDEtools:-dchange(tr, pde, [X, Y, Z, T, U]);

Divide the result by the coefficient of diff(U(X,Y,Z,T), T) and simplify

coeff(lhs(%), diff(U(X, Y, Z, T), T));
simplify(%%/%);

We arrive at:

You want this to be in the same form as the original equation.  You don't want a=1, so we need the exponents to be zero, that is:
eq1 := -Zeta-beta+alpha = 0;
eq2 := -beta+alpha+mu+nu = 0;

This is a system of two equations in five unknowns. You may pick three of the values arbitrarily and calculate the other two.

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.