The GUI does not actually use a Maple command, when it adjusts the plot's characteristics via that right-click menu. What the GUI has at hand is the PLOT call structure.

There is some new-ish "re-draw" stuff related to a command that the Maple plotting commands sometimes inject as text into the plotting structure -- eg. PLOT function call, as structure. It can use that now for some zooming. For example you can zoom out of the following plot, and magically the (post Maple 2022) GUI will cause the plot command to be reformed with a widened range, and reexecuted.

P:=plot(sin(1/x), x=0..1, gridlines=true):
indets(P,"input"=anything); # outputs the following:

{"input" = [table([1 = plot, 2 = [sin(1/x)], 3 = (x = 0 .. 1), 4 = (gridlines = true)]), "originalview" = [0.0000192831942484872343 .. 0.997499999997671694, -1. .. 1.], "originalaxesticks" = AXESTICKS(DEFAULT, DEFAULT, _ATTRIBUTE("source" = "mathdefault"))]}

So, that is actually stored in the PLOT function call (data structure) that the Library plot command produces, in modern Maple. And when the rendered, inline plot is zoomed out then the GUI is now able to send a modified reconstruction of the plotting command call, to the kernel. The original structure only had data computed for x=0..1, but after zoom-out the reexecuted call can recalculate for the wider range.

But the GUI's knowledge of how to adjust the call is limited. It knows how to change the range, but not everything. And not all plotting commands store the details in the structure. And, importantly, there's no place for the GUI to hand back the modified command (or even the new output structure) to the user.

A plot in a PlotComponent does allow the user to access the GUI's stored structure. If you right-click the above plot when inside a PlotComponent, and change the color, then you can programmatically extract the modified plot structure from the embedded component (GetProperty). But unfortunately the spiffy zoom-out redraw/reexecute thing doesn't seem to work in a component. And the GUI lacks knowledge about how to itself adjust the relevant stored command data, in most cases.

That leaves the PlotBuilder. You can issue,
   PlotBuilder(lhs(eq));
and, in the right-panel, modify the range and the gridline choices, etc, and toggle on the "Show Command" button. That can be a handy way to explore various option choices, and see/copy the underlying command. Unfortunately it doesn't yet support discrete data (eg. your list of numeric root values, for which pointplot could be a choice).

@salim-barzani In the attached variant, the lastest parameter values used in the exploration are available (outside the exploration). They are assigned in a list to the name last.

graph_in_simple_shape_accc.mw

You could, of course, subsequently programmatically substitute those values into your M expression, if you also need that.

Alternatively, the following variant has the latest instance of M (substituted) assigned to the name last.

graph_in_simple_shape_acccB.mw

Download r01.mw

You can supply 

    output=L

to the,

   CodeTools:-Usage

command call. The list L can contain items like, value, realtime, bytesalloc, etc.

That makes it return an expression sequence of all those output choices. So you can get your hands on the timing/mem results programatically, as additional returned results.

You can then store items in that returned expression sequence, as you see fit. And your computation's own result is included in that sequence, if value is an item in the list L.

That's in contrast to the default behavior of a merely printed message containing the various stats, and only the value as single returned result.

Is this something like what you're after? (I realize that the final construct is not in terms of functional differential operators.)

I've used PDEtools:-undeclare(prime) to get the derivatives displayed in a "subscripted" notation.

Download BE_ac.mw

I realize that this has some flaws; it can't get by discontinuities in the numeric solving, due to the set-up. But perhaps it might inspire you.

(I can't help but feel that I once knew some better way, but right now cannot recall it. I'm sure I'll kick myself, shortly.)

Download Dipole_fields_ac.mw

Is this the kind of thing you're after?

Download test_sol_for_PDE1_ac.mw

You have mixed up the order of the arguments. You've used a different order, in two different places.

You have DOE set up with the columns (apparently) being:
    epsilon, A, Br, N, beta, m, lambda0, lambda1

But you constructed your procedure Fop to treat its arguments as if the order were,
    A, Br, N, beta, epsilon, lambda[0], lambda[1], m

I'm guessing that the wrong one is your construction of Fop, because there you just used a lexicographic order, and ignored how you set up the columns of DOE.

So I suppose that you wanted Fop more like this: Answer_not_same_ac.mw

Such a placeholder is not required. And you don't need to substitute into a diff(...) result.

Instead, you can just use D(W). You can use it as D(W)(u(t)) in your DE for dsolve. And you can use it for D(W)(u) in your fieldplot call.

And you also don't need to define W the same way, a second time. Since the original W was defined before the parameters were assigned values then calls to W or D(W) will pick up the latest parameter values assigned.

So your same processes can become more tidy and legible, eg, proteins_ac.mw

If you really want to visualize the latest parameter values in the formula then you can just invoke it to print it, eg.
   D(W)(Un);

If you want you can assign D(W) to some name, to avoid having D repeat effort. Eg,
   dW := D(W);
and then use that. Eg. proteins_acc.mw

You could put all your procedures into a package module, and then LibraryTools:-Save that instead.

See also the relevant section of the Programming Guide.

Here are various ways in which one may leverage the collect command to attain the target form.

Download collect_I_ex.mw

You can simply pass the explicit option, to get the result in terms of explicit radicals directly from the solve command.

Download solve_sys_expl.mw

In other words, you don't need to use two different commands to get the roots in radical form.

For a polynomial system involving more than one variable then results involving roots/radicals are not expressed in explicit form by default (even for this low degree).

See also the subsection Controlling the Form of Solutions in the Help page for topic solve,details . (Unfortunately, the difference between a single polynomial in one variable and system is not made entirely clear there, in this respect.)

First, and most importantly, notice that,
   0 < x and x <> 2
evaluates to just,
   0 < x
since <> and = evaluate as if under evalb, with active and.

In this case you are running into the effect of,

evalb( x <> 2 );
             true

So 0 < x and x <> 2 becomes 0 < x and true, which becomes just 0<x . Not what you intended!

This subtlety has caught out many Maple users; you're in good company.

Instead, you can use an inert form such as (captitalized) And .

Download solve_ineq_1.mw

Now, we've obtained more consistent results. And if you're wondering why the slightly different forms of input produce the different forms of output, well, it's because it's more useful to have the choice. (It's been noted previously, and reported, that documentation of all the input/ouput forms from solve are under-documented.)

I am guessing that you'd like all the real roots to be stored, in the case that there are multiple real roots for a given set of parameter values. If not you can pass the option maxsols=1 to fsolve.

Download fsolve_oper_expr_ex.mw

ps. Some people prefer the alternate name MakeFunction instead of unapply. That's moot if you want all the roots, since using expression-form for fsolve is more efficient in that case.

pps. I put the roots in a list, since in the case of multiple roots I expect you'll find it slightly more straightforward to deal with a stored list rather than deal with stored expression sequences.

I hope you realize that numerically integrating the numerically computed derivative is (numerically) generally a bad idea.

Also, I suppose that you realize such numeric computation is not all necessary here, and that (as you say) you are trying to establish something about the three numeric solvers in play here.

But you have been unclear as to why you think that this methodology shows something significant and what precisely that might be. (Also, why differentiate numerically at all?) There are so many controls in play here (eg. the working precision in each of the three numeric solvers in play, tolerances, etc.) that you blanket/vague goal lacks focus and thus meaning. This is the kind of thing where it's more sensible to start with the numeric requirements, explicitly. Then see if you can attain them. Otherwise (IMO) you risk flailing more with the functionality choices and missing something key.

With all of the 1) numeric pde solver, 2) numeric differentiation outside that pde solver (ugh), and 3) numeric integration then any numeric difficulties are going to be a causality muddle. How would you ever decide which might be some weak link? (especially considering that not all the parts are necessary...)

It's not clear to me what you mean by, "...the global form of the conservation law over time".

Is this the kind of numeric integration you were trying to get at? (I'm not suggesting that this is the plot you're after. I include it merely to go with that integration result at the particular t=0 point.)

Download pde_num_hmm.mw

ps. Here's another syntax for the integration and plot above. You original syntax was not ok, using an fdiff call like an expression alongside the (operator style) a..b for the range of integration.

evalf(Int(x->fdiff(sol_proc(x, t), [t = 0]),
          a .. b, digits=10, epsilon=1e-5, method=_Dexp));
plot(x->fdiff(sol_proc(x, t), [t = 0]), a .. b, size=[500,200]);