@vv Yes, I see that my dsolve solutions eventually (as the parameter increases) are no longer roots of the polynomial. Thanks for pointing out my error. Yet the parameter k is a coefficient of the polynomial, so the roots are continuous (complex-valued) functions of k (right?). How can those functions be plotted?

@mmcdara I agree. What Maple calls Quantile(X, ...) is commonly referred to as the "inverse of the CDF of X". Isn't that correct? I'm sure that this is what the OP meant; nonetheless, you're right to point out that for many distributions X, "inverse X" means something completely different from "inverse of the CDF of X".

@Racine65 I'd begin with a photo of a teapot -- an orthogonal projection of a profile view. Then trace or freehand the projection only (in other words, a 2D image -- a silhouette or shadow) onto graph paper. Decide what features you want to represent -- let's say handle, main body, lid, spout. We'll make separate Bezier curves for each. Use the graph paper to get numeric values for the control points. Note that the two endpoints of a Bezier curve are usually the only control points actually on the curve. Don't worry about the number of points per curve. Make tubeplots of each. For the body and lid, the Bezier curves will parameterize the radius to produce surfaces of revolution. For the handle and spout, they will parameterize the central axis. Then put the plots together with plots:-display. Then dress it up with standard 3D plot options for coloring, glossiness, removing the gridlines, etc.

@Carl Love I think that if there's going to be a finite scale of zoom factors, then it should be a geometric scale (with a ratio about 1.1) rather than the current arithmetic scale (with a difference of 25%).

@dharr Here is a better image

I do (in Maple GUI) "Right click on plot border" --> Copy Special --> Copy as Image, then (in MaplePrimes) Ctrl-V.

@Carl Love In an earlier version of the code above, I forgot to include the binomial coefficients. The code above has been corrected.

@Racine65 If you want to get a teapot, consider using plots:-tubeplot rather than plots:-spacecurve. The radius of the tube can be set or parameterized.

@Chouette There's also evala(Normal(...)) and radnormal(...).

@acer It should be noted that parametric returns nested lists, which need to be flattened for most programmatic uses, so they're syntactically different from other solve return values. It's not a big deal here. However, the more-usual case is that you get a piecewise of nested lists, or even a nested piecewise of nested lists. They are awkward to deconstruct.

@ecterrab Could you please show how to apply it in this case, where the goal is to reduce the independent variables to simply y? I see that ?dchange has an example of reducing the number of independent variables, but I can't get anything to work in this case (which seems like a simple case, since my one line 

convert(simplify(eval(PDE2, u(t,x)= U(x-a*t)), {x-a*t= y}), diff)

does it). It would be preferable to have one command, dchange, to do it rather than my composition of convert, simplify, and eval.

@acer That the first workaround that you proposed -- adding a dummy variable to the solved-for set -- works provides a small bit of indirect evidence in favor of my theory. 

@Kitonum I think that a slight improvement would be continuity with respect to time. You can do this by making the t-range an integer multiple of the period.

@sand15 There's no need to normalize your random numbers to a 0..1 float scale. The integers generated by naked rand() are good enough, and more efficient.

@Magma In most cases, randomize() should only be used as a top-level command, and it should only be used once (per restart).

Kitonum's code is extremely inefficient. I will try to rewrite it.

@Magma Your results indicate that the output= permutation part isn't working. I'll come up with something else in a few minutes.