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]);

Are you perhaps looking for something like the following:

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

I'm not 100% sure how your code is supposed to work.

On my machine the following takes about 14sec, compared to the original taking about 24sec.

Formula_II_acK.mw

(It's not easy to work with this code when it's "behind" a Button. But I've put my edits back there because that's how you gave it to us.)

A significant portion of the outstanding time is now taken by the coeftayl calls.

I hazard a guess that some of the following might be possible:
a) Form parts or a variant of `Eq` just once, and then instantiate more quickly when the input Components's values change.
b) Figure out a way to have the coeftayl results be done symbolically, so as to instantiate a template rather than do its heavy lifting for each set of Component values.

I expect that there are further improvements possible. Alas I'm out of time for this, sorry.

I noticed that if I changed parameter N from 0.5 to 1/2 then the solve step seemed to slow down. Hopefully that is all "above board" and not any hint that your current results are actually slightly suspect.

Here I use the formula as in your picture, with the x- and y-components switched. And so a similar plot can attain.

plots:-display(
  plots:-spacecurve([t*cos(t),t*sin(t),2*sqrt(2)/3*t^(3/2)],t=0..Pi,
                    thickness=4,color=black),
  axes=normal,view=[-3..3,-2..4,-2..8], scaling=constrained,
  orientation=[50,70,0]);

You also mentioned 2D plotting, but that could mean several different things and you haven't explained what you are after. Perhaps a projection of that spacecurve, onto the x-y or some other plane?

Download poly_puzz.mw

The roots are R, and the polynomial is P.

@C_R Notice that your code does not assign the expression to a name, and merely reuse it at the different values of Digits.

Rather, your code explicitly re-forms the expression at the different values of Digits.

The 1e-11 invocation forms the very same float expression at both Digits settings. But the 0.1*10^11 forms different float expressions at those two different values of Digits.

The differing coefficients makes the expressions being passed be to is be different (ie. different address, etc).

That makes the remember table get used for one flavor (where the expression passed is identical) but not the other (where the expression passed is different).

Download flt_ex.mw

On the other hand, if you formed the two expressions just once, up front, then the remembered can interfere (here, consistently).

eval_of_range_in_ln_ac2.mw

So, the cause is that (only) one of those constructions produce identical results at differing Digits.