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. 

In general this would just be a matter of adjusting the view option, giving both a horizontal and vertical range.

But I think here you want to see the intersection point between a line through the points [30000,0.2838019625*10^8], [60000,0.2122282820*10^8] and a line through the points [30000,0.2838047526*10^8], [60000,0.2126554523*10^8].

The y values are identical for the first few decimal places, so it will not be possible to visually dstinguish them.

In your first attempt you pass an expression to plot3d with the inverse fourier transform not evaluated, so at each point in the plot there will be an attempt to find it, which will fail. The solution is just to evaluate it before the plot, which can be done by assuming t>0. For inverse transforms you almost always need assumptions. (I added with(inttrans).)

ft1.mw

I don't quite understand point 1. Perhaps you want to be able to choose other values of S__1 and S__2? I would set

params := {S__1 = -1, S__2 = -1};

then use eval for the things you want to use these parameters in, e.g.,

eval(P, params);

This way, you only need to change the params line to run the worksheet with your new parameters.

For point 2, I'm not sure if you are asking how to find the conserved quantity, which I've shown how to do in the attached. If you are asking about your empty plot, that is because you need wider contour intervals, which I adjusted. You may want to try without the "contours =" option to see which contours Maple comes up with first, then if necessary make a finer adjustment by adding the contours option.

L1C.mw

It would be helpful if you used the same notation as the paper. Here's my attempt, but I don't get the result in the paper.

Download C1.mw

In this case it's probably easier to do the phase portrait with DEplot.

Download conserved.mw

dchange from the PDEtools package can transform from one set of variables to another.

Download polar.mw

I agree it is a bug; it is still present in Maple 2025.1.

Edit: SCR submitted.

A simple workaround is just to substitute a non-indexed name before carrying out the command, and then reverse the substitution on the result; if you have more complicated requirements, I am sure there can be a more sophisticated solution.

Download PolynomialToPDE.mw

When you wrote puzzle(n) := etc you defined a procedure - notice that the output says puzzle := n -> etc. But you did not invoke that procedure - if you do you immediately get the answer. (I would define  puzzle by puzzle := n -> etc, since then the defining syntax is distinct from the calling syntax, but modern Maple allows both.)

Download test.mw

@kencom1 In Maple 2025, your worksheet does not run as you show it; the notation f[0](1) when setting the initial conditions leads to an error. You are having problems using the same name for procedures and indexed names, which become tables and overwrite the procedures when you assign to them. The solution is to have different names, e.g. f for your procedure and F for your initial conditions and odes;

Download notes.mw