@student_md Here's a module to create your matrix:

MyOuterProd:= module()
local
   phi, t, #function and argument
   K,
   a:= j-> doublefactorial(2*j-1)/j!,
   g:= (m,s,j)-> a(m-j)*a(j)*a(s-j)/a(m+s-j)*(2*m+2*s-4*j+1)/(2*m+2*s-2*j+1),
   multiply:= proc(n,m,l,s)
   local j; 
      if n <> l then 0
      elif m > s then thisproc(n,s,l,m)
      else 
         sqrt((2*m+1)*(2*s+1))*2^(K/2-1)
            *add(g(m,s,j)/sqrt(m+s-2*j+1/2)*phi[n, m+s-2*j](t), j=0..m)
      fi
   end proc,
   ModuleApply:= proc(M::nonnegint, K::posint, phi::function)
   local M1:= M+1;
      thismodule:-K:= K;
      thismodule:-phi:= op(0,phi);
      thismodule:-t:= op(phi);      
      Matrix(
         (M1*2^(K-1))$2, 
         proc(i,j)
            local n,m,l,s;
            n:= 1+iquo(i, M1, 'm');  l:= 1+iquo(j, M1, 's');
            multiply(n,m,l,s)
         end proc,
         shape= symmetric
      )
   end proc
;
end module
:  

The parameters are M, K, and phi(t); for example

MyOuterProd(2, 2, phi(t));

Please check the output carefully for several examples to make sure that it's what you want.

@vanzzy In cases where there are no possible explicit symbolic expressions for the roots of a polynomial (as is usually the case when the degree is greater than 4) Maple uses expressions RootOf(...) where each of these represents all of the roots of the polynomial in its parentheses. The command allvalues will split these RootOfs into indexed expressions (usually of the form RootOf(..., index= n)), each of which represents a single root of the polynomial. If you then supply values of the parameters (such as your k) and use evalf, you'll get decimal approximations for all the roots. So, using your example above, do

H:= allvalues(H);
nops(H);

and you should get 5.

Note that my Answer below that uses dsolve to plot the roots is incorrect.

@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.