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

Here's how I would follow the Youtube simplifications in Maple land. I also show that term1 and term2 are the some parameters. Edit: I mistakenly thought that @Carl Love was using t for x, but I now see he was referring to a different post.
 

Download MapleLand.mw

OK, I now understand you want the plot of u(y) to look like the one you just posted, with a slope of zero at y=0 and going down to u(1) = 1.
Now psi = D(u), but in the odeplot that you labelled u(x,y) you asked to plot D(psi), which is not u. u is instead the integral of psi.

pltN[j] := odeplot(
        dsol[j], [y, D(psi)(y)], y = 0..1,
        labels = [typeset(y), typeset(u(x, y))], labelfont = ["TIMES", "italic", 18],

So just add another ODE to your system

diff(u(y),y) = psi(y), u(1)=1 

and plot u(y).

Edit: Of course you will have to change your boundary conditions that assumed D(psi) = u instead of psi = D(u)

Yes, that seems like a problem. Using LinearSolve is faster. (If you can use a float calculation NullSpace is much faster.)

Download NullSpace_vs_kernel.mw

Edit: If you want it in the same format as NullSpace:

NullSpace.mw

What do you want to plot? You could use plot3d to plot abs(FF) vs x and P1 or Re(FF) vs x and P1 or Im(FF) vs x and P1 or argument(FF) vs x and P1. Here I use complexplot 3d to plot abs(FF) coloured by the phase. The plot is not very informative, probably because your function has numbers with wildly different exponents, e.g. 10^(-2750) and 10^5 but these numbers have only about 5 digit accuracy. Probably you should go back to a version in symbolic form and then evaluate to numeric form with a high setting of Digits so you don't lose accuracy.

PLOT11.mw

Edit: If you really want to find where the function is zero with implicit plot, you should explore more limited regions - the plots above can be a guide for this. I suspect again the quality of your function is making it hard for the internal solver.