Yes this factorial in the denominator jumping drives me crazy. It seems to be a side effect of the new 2-D entry system that abbreviates entering x^2 cursor-right + 2 to  x^2+2. It is still present in Maple 2025.1, so there is currently not a solution.

Well, this is probably cheating since I saw the answer, but maybe I would have thought of it anyway.

Download test2.mw

Of course it is a bug, but notice that solve gives the same result. It's the classic problem of introducing extraneous solutions by squaring both sides of an equation.

Download realsolve.mw

[Edit - modified so step 1 works as described in the literature - see more examples with refs below.]

Download AllstepsBurgerNew.mw

Example 2: OP's Hirota Bilinear eqn

AllstepsBurgerNew.mw

Example 3: KdV equation

AllStepsKdVnew2.mw

If you create the Matrix with option shape = symmetric, then you will get only real eigenvalues and eigenvectors, and Maple will use a more efficient method specialized for symmetric matrices.

(If you just want to get rid of the +0*I after the fact, you can use simplify(..., zero)).

You didn't say how you want to export your plot, but if I assume you want vector graphics then there are two possibilities, encapsulated postscript (EPS) or SVG. (possibly also PDF but I've never tried it)

If I make a plot with no background, say plot(sin(x),axes=none): and export it as EPS with the right-click menu, then look at it with CorelDraw I find it has a white background (which I usually delete and then process further), so this does not work for you.

However, exporting as SVG gives a file that does not have a background, just the curve, so this is what you want.

If you actually want to export as a bitmap image (with loss of resolution), then I don't have much experience with that. Which formats would you want then?

An extract from the ?updates,Maple2017,Performance help page: "Maple 2017 includes a new compiled C implementation of the FGLM algorithm for computing a lexicographic Groebner basis from a total degree basis when there are a finite number of solutions.  This routine is used automatically by Groebner:-Basis when applicable and by the solve command when solving polynomial systems.  The example below runs about 200 times faster in Maple 2017."

The ?updates,Maple2018,index help page and the current help page for ?Groebner,Basis say "The Groebner[Basis] command was updated in Maple 2018", which I assume means the command itself has changed but not necessarily the algorithms.

(If you are using the bases for solving systems of polynomial equations, then there are new methods that don't use Groebner bases that may be more efficient.)

If I understand correctly what you mean, use ctrl-shift-K to insert above and ctrl-shift-J to insert below. (command-shift-K and command-shift-J on Mac). See ?worksheet,reference,hotmac or ..hotwin or .. hotunix.

The default behavior was changed in Maple 2020, see ?updates,Maple2020,IntegralTransforms.

To get the behavior you want, use setup(computederivatives = false);

Download laplace.mw

This is an answer to the followup question, now deleted, about plotting the solution in cartesian coordinates. Here is an example:

PDE_upd_5.mw

Note that VectorCalculus:-Laplacian(u(xi, eta), 'bipolar'[xi, eta]) assumes a=1, so I substituted your earlier laplacian.

Since you have already the functions x(xi,eta) and y(xi,eta), you can use the parametric version of contourplot.

PDE_upd_5.mw

I would probably play around with custom contour values because your large spike dominates and so the contours near zero do not have good resolution, but the I think you cannot then use the colorbar.

NOTE: I replaced VectorCalculus:-Laplacian(u(xi, eta), 'bipolar'[xi, eta]) with your earlier hand-calculated Laplacian, since the built-in one uses a = 1, which you had before but is not applicable here. VectorCalculus is supposed to have a mechanism for changing the a parameter (SetCoordinateParameters) but I couldn't get it to work with Laplacian. 

Yes, your substitutions missed an x. You can capture this with 

eqs := {coeffs(collect(numer(normal(eq)), {X, Y, Z, x}, 'distributed'), {X, Y, Z, x})}

but I think that you cannot constrain tanh(k*x + p*y - t*w) together as one variable, since x, y, and t are independent. If the ln was not there you could expand(tanh(k*x + p*y - t*w)) and substitute the individual sinh, cosh etc, or just convert to exp, but in both cases the ln is a problem.

I would almost always use simplify for this sort of thing, not is.

Download distances.mw

There were errors in V. Now it works but is very slow, even with only a 10 x 10 grid. There are common expressions in the integrand that some hand coding could optimize. One could also play with some of the integration options (method, errors).

Download fokas_method.mw

Here's an attempt at a ternary plot. If you want tick marks along the edges of the triangle, that would be tricky to do, but undoubtedly possible. (One edge would be easy because it is along the x-axis.)

Download trinomial.mw