@Khair Muhammad Saraz 

Please attach your worksheet.

Taking nops of a Matrix or a Vector (one of its rows or columns) does not return the number of elements. Your MM is a Matrix. The command for that is numelems, not nops. There are better ways to get the dimensions, IMO.

Perhaps you intended on of these:

listcontplot(MM[..,2..], axes=boxed);

listcontplot(MM[..,2..]^%T, axes=boxed);

listcontplot([seq([seq(MM[i,j],
                       j=2..op([1,2],MM))],
                  i=1..op([1,1],MM))], axes=boxed);

Let me know whether you'd rather do interpolation (linear, spline, cubic) instead, to get smoother curves or nicer axes tickmarks.

@MaPal93 

Regarding your attachment above that balked after 55 iterations (I think maybe it was named fsolve_does_not_work.mw):

Once INDEX gets up to about 0.55 (or 0.554...) some of the solved variables become very large (in a solution). Then you have both small and large values in the "root". In order for fsolve (or you) to establish that it really has obtained a root it performs a residual calculation. But as things get worse the required working precision (Digits) needed to do that accurately gets larger too.

Experiment indicated to me that default Digits=10 was producing bogus results already. Still indicative that solution values might be approaching a singularity. But not accurately. Raising Digits (30, 50, 100...) only brings temporary relief.

Moreover (and I know this sounds facile, but I'm being serious), it can help to pause and try to define the notion of a "root" -- in this float context. Once you have such, you can program for it.

By default the residual check (after resubstitution; forward error) looks to inbound Digits to fsolve for its scale. You can make an end-run around that by using operator-form instead of expressions, for the equations. Such black-box procedures are free to temporarily raise just their own working precision, for just each time they are invoked. This is an end-run around the catch-22 that raising Digits at the higher level -- as a means to evaluate the expressions more carefully at each iteration -- also has the effect of making the residual checks more strict. Using black-box procedures can separate the working precision for the functional evaluations (Digits now set withing the scope of the procedures) from the accuracy requirement in the residual check.

Which of the names I1,I1,alpha,beta,y,d,h,Q1,Q2 are the parameters over which you want to minimize, and which have (some of the) numeric values that you've used?

Note that the objective TR you've used was defined using y = 0.7, d = 2*10^9, h = 35, k = 10^7 and alpha = 0.02*exp(0.4*I1), beta = 0.02*(1 - exp(0.5*I2)).

So why then do you include constraints that only include alpha,beta,y,d,h ?

You also defined TC and TC2 using Q1 = (alpha + beta)*d, Q2 = x*d*h - delta*(alpha + beta)*d*h.

So it's quite unclear what role you intend for your constraints, and what the independent parameters are.

@MaPal93 I don't know that it could be done using the usual substructure generated by the legend item.

Of course, it would also be possible to fake the legend by putting three straight lines in a polygon that is placed carefully in an otherwise blank region of the plots interior. But that kind of ad hoc solution requires effort (locate such a region, make a nice scale, etc.) The up-side is that you can make such lines as long as you'd like within the restrictions of the region.

@Khair Muhammad Saraz To address two more followup queries by you:

For a 2D contour-plot from the Matrix MM you have two choices:
1) Do it directly, using the plots:-listcontplot command.
   (That's easy to do directly, but the axes tickmarks are trickier
    to get to match some arbitrary ranges). Eg,

   plots:-listcontplot(MM[..,2..],tickmarks=[[],[]]);

2) Using Interpolate, and then use plots:-contourplot on that.
  (The first part is trickier, but has also been done before by me on
    this site. There's also an Example worksheet on Interpolation
    and smoothing that illustrates the first part.
    See also,
      Help Topic, ?examples,Interpolation_and_Smoothing
      some old post)

In my Answer's original attachment I already showed you how you could make nice tabulated data from a Matrix like MM.  dq_code_1_acc.mw

@nmacsai The baseline for the color-bar doesn't appear to adjust properly w.r.t. the bottom of the title (in the way that the 3D portion does, say).

It seems that enlarging the view of the COLBAR's contents, along with judicious use of size for both parts, can make a suitable vertical adjustment.

Can_I_change_the_location_of_the_color_bar_caption_in_Maple_2024_sol1_ac1.mw

Unfortunately not all of that was automatically determined. Automatic determination of the adjustment values might be possible, using font height (if obtainable in pixels, which can be an OS thing) and the view-to-plot-region-size ratio.

ps. Combining by plots:-display may well strip off such size adjustments. That's why I did a substitution for COLBAR's elements only after plots:-display.

@sursumCorda On my machine I need all the following for such "subexpression labeling" to occur,

1) interface(typesetting=standard):  # not default for my Maple 2024.0 on Linux
2) interface(labeling=true):  # default
3) interface(prettyprint=2):  # 3 by default

After that, it works. Unfortunately line-breaking of output is disabled by it, and the Java GUI is as good as the old Classic GUI at figuring out a good split. Manually doing,
  interface(labelwidth=<someposint>):
is possible but I find that often no single value works as well as I'd hope.

Your code has a couple of instances of,

   else then

but that's not valid syntax.

Perhaps you just intended,

  else

@Carl Love I wonder whether it could be even more visually informative if the orthogonal trajectory lines were shaded by steepness.

That's easier to do with 3D space-curves. For 2D curves it's less directly easy, though of course there are a variety of ways, eg. transform from 3D, manual construction, etc. (I've come across the need for easier 2D curve coloring before... I should collect prior notes and report suggestions.) 

@C_R 

[edit] For this followup snippet of yours one can examine the code and see that, yes, the name unknown is being used as a local.

    showstat(`type/satisfies`);

Perhaps a developer wanted to make a play on the idea that the predicate foo in a type(expr,satisfies(foo)) call is often an anonymous procedure.

I forgot to show in brief form the effect of CollapseNested. So, on the chance that it might not have been crystal clear to all readers:

The two forms print slightly differently.

The key thing is that the collapsed form is just a single int or Int call, which controls the whole process and can utilize assumptions on the inner variables of integration (if it possible for it to infer it from the nested bounds).

Download IT_EM_ex.mw

I edited my previous Reply to use,
    x::RealRange(Open(Pi/3), Open(Pi))
instead of,
    RealRange(Open(Pi/3), Open(Pi))
in order to show that the only assumption placed & required there is on x. No manual assumption was placed on the variables of integration. (That was already occurring, as explained a few Replies back. But I wanted to show the aspect.)

@sursumCorda Your first example takes quite a while, and emits Warnings.

As I mentioned, using a "collapsed" form may allow a single int or Int to represent the whole triple-integral and, when computing, to pass the outer integrals' range inferences on to the inner ones. Using that form gets rid of the Warnings for your first example.

Using the collapsed form can sometimes speed things up. And sometimes it resolves an example that otherwise wouldn't succeed.

I asked my oracle and the sibyl responded with a simplification chain than can get the zero for your first example.

Download sursumCorda_int_ac.mw

In contrast here's that original first example. Granted, some timing differences may well be due to different form of the integrand (including combined trig, etc).

I found that direct timing comparisons within the same session were problematic, due to remembered results (and intermediate internal results). I found that example even more prone to such difficulty than usual; I don't know if that's due to the nature of the example or something to do with Maple 2024.

@nm The Physics updates package doesn't get installed by default. One has to choose to install it.

@Khair Muhammad Saraz I made two surfaces in each of my latest 3D plots.

Contrary to your apparent new claim, I did in fact use the approximate data for one of the surfaces, each time.

I stored the approximate data in Matrix MM. I used the plots:-surfdata command to render a surface using the data in Matrix MM. (In Matlab you showed using surf. That's similar to using surfdata in Maple.)

Can you see that my latest worksheets have a red exact surface as well as a blue approximate surface? Earlier you wrote that you wanted both. Can you see that the blue surface is constructed using the computed approximate data from your loop?

If you don't want them both together in the same plot then simply don't bundle them together with last display call. They're even already assigned to names which make it clear, Pexact and Papproximate.

This is my last response on this. You really should study the code and try to understand it and learn to program in Maple on your own. This is not unreasonable stuff. Try to understand properly what each part does.

ps. I think that storing the data in a Matrix (like MM) here is easier than storing it in tables (with that awkward copy each time through the loop). Using a Matrix makes it easier to plot afterwards, and (I expect) more understandable if you were to study it.

@nm Is this something (even for only some x...),

-sqrt(2)*EllipticF(sqrt(2)*sqrt(sin(x)/(1 + sin(x))), sqrt(2)/2)