@sand15 That seems correct. I wasn't aware of that option. Is it on a help page? 

That value seems to be unchanged by calls to plots:-odeplot; it must be accessed after calling sol(t) for numeric t.

Almost every time that I've used assign to assign the values of pure symbols (like x, y, z) rather than table entries, I've gone back and rewritten the code to use eval. Assigning values to simple names seems like an easy way to cause bugs when you inadvertently use the same names later expecting them to be symbols.

restart:
gm := V -> 1/sqrt(1-V^2):
T := (1/sqrt(1-V^2))*(w-k*V):
S := (1/sqrt(1-V^2))*(w*V-k):
f := -135/4*T^5+369/16*T^3*S^2+47/4*I*T^4-93/16*I*T^2*S^2+T^3-2/3*T^3*S*B-27/16*S^4*T+3/16*I*S^4-1/3*T*S^2+2/9*S^3*T*B:
Hgen:= simplify(rationalize(f)):

I don't think that setting _EnvExplicit to true can possibly do anything beneficial for a polynomial 
of this complexity. While I also don't think that it can cause harm, since there's no reason to set 
it, I'd rather not risk it. So, the next line is commented out.

(* _EnvExplicit := true; *)
H:=[solve(numer(Hgen),w)]:

If the solve goes as planned, then H will be a list with a single entry, a RootOf expression.

nops(H), type(H, list(RootOf));
                            1, true
So, H is indeed of the the form that I hoped.

The part of this worksheet that solves your issue of the slowness of allvalues is in the next 
line. The key point is that numeric values are supplied for B and V before the expression is 
passed to allvalues.

Hav:= CodeTools:-Usage([allvalues(eval(H[], [B= 1, V= 0.8]))]):
memory used=157.39MiB, alloc change=8.00MiB, cpu time=1.09s, real time=1.12s, gc time=78.12ms
nops(Hav), type(Hav, list(RootOf));
                            5, true

There's no need for you to re-enter all those plot options. Your entire plot can be done as below.
 
plot(
   Im~(Hav), k= 0..1, discont,
   color= [red, black, green, blue, yellow],
   linestyle= solid, thickness= 2,
   axesfont = ["TimesNewRoman", 18], 
   labelfont = ["TimesNewRoman", 22], labels = [conjugate(k), Im(conjugate(omega))],
   axes = "boxed",
   gridlines= false
); 

Download ImagRootPlot.mw

@mmcdara You mentioned the need for a consistent order to allvalues solutions. This is a very complicated issue, and I often find that there's no satisfactory solution for whatever problem I'm working on. First, note that allvalues only has an affect on RootOf expressions. Then very carefully read the help pages ?allvalues, ?RootOf, and ?RootOf,indexed for details on how you can set the order and how you can change it.

As you've been told before, please upload a Maple worksheet showing the error using the green uparrow on the toolbar of the MaplePrimes editor. If the upload gives you an error message, just ignore it: It likely only means that the worksheet can't be displayed in MaplePrimes; we will likely still be able to download your worksheet.

I guess that RK-4 is 4th-order Runge-Kutta. Is that right? What are NSFD and LWM?

@cinderella 

Do

fsolve(Val, Q= 10..20)

to restrict the solutions to that range.

@vanzzy When I said "using your example above", I meant using H exactly as it had been defined in your previous code. That is, you must first do

H:= [solve(numer(Hgen), w)];,

as before, and only then do

H:= allvalues(H);.

The reason is that the allvalues only affects RootOf expressions. If you don't first use solve, then there are no RootOfs.

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