@mmcdara Thanks, that's the kind of thing I'd surmised. But the numeric requirements are imprecise.

My concern is that such a process is not a test of only the numeric pde solver. Its results rely upon not only behavior of the three different & separate numeric solvers I mentioned but also upon their interplay.

Failure to attain some (as yet unstated) degree of numeric accuracy cannot therefore be a proper gauge of numeric performance of any of sole component part of that process. Eg. failure might say nothing about the numeric pde solver alone.

I supposed that accuracy of the integration step would need be coarser than that of the diffentiation step, that in turn being coarser than that of the pde scheme. But the particulars of that might need be problem specific. And attaining an upfront, explicit, accuracy goal might need per-problem tuning at any/all of those levels.

@GFY By the way, it is also possible to do it with fsolve (as a modification/correction of the earlier attempt). Here, I augment with undefined, for the absent root values.

You can also scale any you choose, eg. by 3.

I use option remember to avoid calling fsolve multiple times for the same input (ie. same omega__v, or r, value).

You could make this more sophisticated.

Download fuxian_ac2.mw

@GFY Firstly, note that there is a range/set of values of omega_v for which there is only one positive root Omega of det_num. That is why there is a middle section with no upper root. It's not that the method fails to properly compute them; they don't exist in that region.

[seq(RootOf(det_num, Omega, index = i), i = 2 .. 3)];
plot(%, omega__v = 0.95 .. 1.01, size = [400, 400]);

Secondly, you are not being clear about what you mean by the term "root". Even if you were able to construct explicit formulas (with radicals) those need not match up with the apparent "curves" in the way you expect. (That can happen even for a parametrized quadratic.) The numeric rootfinding method in play here involves generation of all the roots (per omega__v value), which you may sort. That's what you had originally; you computed all roots (for each given omega__v) and then you applied sort. And that's essentially what evalf/RootOf does here too.

It's a mistake to think that I did something fundamentally different than what you were trying earlier. Using indexed RootOf was just easier than repairing your procedure to inject NULL in the cases/region of less actual positive roots.

You might be able to do numeric rootfinding by a homotopy (or other) approach instead, which gets curves of solid color. If you want color to be preserved across some crossing point of roots then you should at least define that concept.

Another way, with uniform color,

plots:-implicitplot(det_num,omega__v=0.0..2.0,Omega=0..2,
                    rangeasview, gridrefine=4, signchange=false);

You might be able track a "root" curve through a simple crossing/switch by matching derivatives.

@maple2015 You can alter that line to,

   newintg :=  collect(intg,lhs~(inter),uu->simplify(uu,size))

so that lhs~(inter) becomes the list [r,t] or [u,v], etc, as appropriate.

ps. In your original worksheet case 3 had a different substitution, and the numeric integration did not complete successfully on my machine. If you have a revision, you could attach the actual worksheet here.

I have deleted a duplicate Question thread on this.

If you have followup queries on this, then please put them here instead of spawning wholly separate new Question threads.

@Carl Love It's available by doing right-click on his RUN button and choosing the item "Edit Click Action".

Part of the time is spent in a numeric integration in variables r & t. It seems faster to do,
  newintg := simplify(intg,size):
nd then,
  evalf(Int(newintg,inter,digits=15,epsilon=1e-7,method=_cuhre))
without setting Digits=7 (it's Digits=30 earlier, though I don't know exactly why).

It may also be quicker to do,
   newintg := collect(intg,[r,t],uu->simplify(uu,size)):
than to hit the whole with simplify(...,size).

The other part of the cost (about half, maybe) seems to be initial construction of formulas, in particular the line,
   st := [seq(coeftayl(Eq, y = 0, j), j = 0 .. M-2)]
but I haven't looked at that.

@Susana30 You still haven't explained what 2D plot you want.

Maybe it's one of these...
some_plots.mw

(check for mistakes. I don't really know what ranges you want, or if any restricted.)

There are several different ways to construct these (or project them). I duplicated some computation, to illustrate.

@Andiguys 

Your worksheet has s and s_1 as tables. And you know the tau0 values used. (You can also access the indices of the tables programmatically.) So constructing either a mx2 Matrix or a list-of-lists should be straightforward, and you can pass such into the plot command.

Eg,

plot([seq([uu, s[uu][1]], uu = sort([indices(s, nolist)]))])

In my Maple 2019.2 the optimal values for s_1 don't really vary (except for some cruft artefacts of the HFloats displayed). I doubt that it's working as you intend. Increasing Digits:=30 at the very start doesn't help much. It's not clear to me whether tau0 is actually playing the actual role in your objectives that you want.

I'd suggest that the optimization itself might need fixing here. But it lacks comments and proper explanation.

   for tau0 in {0.3, 0.35, 0.4, 0.45, 0.5, 0.6} do
     ...
   end do

@Ronan I'm glad if things are working for you.

Consider the following example. The global name typeset is not a protected name, so you can assign it a value. But if you did so then you'd also need to uneval-quote it when utilizing it in that plotting way, so as to not accidentally get its value.

restart;
plot(x,caption=typeset("hey, ",x),size=[400,200]);
typeset:=33:
plot(x,caption=typeset("hey, ",x),size=[400,200]);
plot(x,caption=':-typeset'("hey, ",x),size=[400,200]);

So, the fix is to use it as ':-typeset', similarly to what's done in your code for several other unprotected names used for options, etc. That should appease maplemint.

The phrasing of that mint warning is slightly confusing. Most often (in practice...) it appears in cases just like here (ie. your :-typeset usage) where some name is wanted as just the name itself, not any assigned value. The other case -- that's less common -- is where one intends to utilize it as a global variable but didn't declare is a a proc global; pretty nick-picky if not being assigned, IMO. So the warning is slightly misleading. Indeed it seems to have led you to guess that the fix in this situation would be to declare it as a proc global. But that wouldn't be right, and wouldn't help in the problematic case of its being assigned a value.

In summary, if you see that warning from mint/maplemint, and didn't intend to use the name for its value, then check to see whether you ought to quote the global instance(s) in the code. Trying to break your own code -- here, by assigning some value to such names -- can be a worthwhile exercise.

Note the effect of forget here. (The `is/internal` procedure has a remember table, for example.)

Download is_remember.mw

The is command can leverage some greater working precision, internally. That can be due to its trying to deal with numeric issues (eg. roundoff, loss-of-precision, etc) at the inbound Digits setting. It doesn't always suffice.

@segfault No, you don't have to do hand-edits. Doing it programmatically is trivial.

I made no suggestion such as "...to take the limit of the function f[4] itself". That's nonsense. Member Hullzie did not suggest that either.

Since f[4](0) is not a name in the technical Maple sense, you can temporarily replace it in your expression. I referenced three programmatic ways to accomplish that.

Hullzie suggested replacement by the name f[4], but that suggestion has quite obviously confused you (perhaps due to the re-use of the indexed name). But that suggestion would not compute a limit of f[4] itself, despite your claim. And you're also free substitute with some other name. ...Or freeze that problematic function call, or frontend the limit.

What Hullzie showed was that if you replace the unevaluated function call f[4](0) by an actual Maple name then the limit you want will attain. You don't have to rewrite the expression or hand-edit it, to effect such a replacement.

The Help-page for limit clearly shows that a name is expected for the left-hand side of the second argument.

f[4](0) is not a name in Maple.

There are many similar instances in Maple where names are expected on the lhs of procedure parameters. Eg, this won't work either:
   plot(f[4](0)^2, f[4](0)=-1..2);

One can freeze the term in question, or substitute by an actual name (symbol, or indexed symbol) of one's choice, or frontend the operation.

@Susana30 An alternative,

Download shb2.mw

@minhthien2016 Using that same procedure F from my last response above,

Download minhthien_ts2b.mw

ps. Since you haven't provided an explicit characterization of the goal then I'm finished here.