@tomleslie This is a correct answer.

@mmcdara This a relatively longwinded and inefficient way to accomplish a simple task.

I'm sorry, but it should be mentioned, because users should not be encouraged to utilize set results (eg. from fsolve or solve, etc) in this way.

@AHSAN You don't seem to understand. Your explanation is incomplete.

What do the colors and gridlines mean, in your latest image?

How do the few fixed pairs of values for two parameters beta and lambda figure into anything like a image with a continuous gradation of color?

There is no point in merely repeating the same inadequate explanation and image. That's not very respectful. You have to explain in detail the precise connection between your mutivariate expression and the plot you want.

I am finished here. 

@AHSAN I am waiting for you to provide a proper description of your goal.

As far as I can see your latest description and image, have little (if anything) to do with your original Question and  attachment (which needed only a few corrections).

I think that guessing what exactly you might want is unproductive.

@J F Ogilvie For what magnitude of relative accuracy are you looking?

@AHSAN Your question asked about a 3D plot, with some like a legend inside some blank space (eg. within the axes).

This image you just added above is not like what you asked, at all.

You should state clearly and precisely what you want  when you first ask your Questions.

@tomleslie The OP's equation was given as,

   64/19683*t^9 - 16/729*t^7 = (-t + sqrt(t^2 - 4))/2

which is equivalent to your expr2 (though you happen to have the radical isolated on the rhs). You squared both sides of your expr2, and expanded, to obtain expr3. But not all real roots of expr3 are actually roots of expr2. In particular the approximate root 2.1213203436 of your expr3 (or your expr4) is not a root of your expr2 or of the OP's equation.

The only real root of the OP's equation that is greater than 2 is given by the approximation 2.5579901741.

Also, if only real roots are wanted then there's little point to supplying the complex option to fsolve, only to have to further select the purely real roots from those.

I understand the appeal of creating your expr4, since as a univariate polynomial fsolve will find all roots for it, automatically. But the roots of your expr4 would need to be each checked for validity as roots of your expr2.

If the OP's equation is assigned to eqn then the fact that it has only a single real root greater than 2 can also be gauged by,

fsolve(eqn, t = 2 .. infinity);
fsolve(eqn, t = 2 .. infinity, avoid = {t = %});

@grheppler The first link in the Answer above starts off with this bullet-point for OSX,

I put the emphasis on that word Preferences, to distinguish it from the word Options.

@TechnicalSupport You're suggesting the calling sequence for plotting two scalar expressions (resulting in two curves). It seems clear that that is not at all what the OP means.

The OP has used the syntax for plotting a single curve, using the parametric calling sequence. It seems clear that the OP understands this and is choosing this deliberately.

That calling parametric sequence does not (yet, Maple 2021.1) support units. That's the OP's problem and complaint.

@Jean-Claude Arbaut This multiplication crash for member vv likely has nothing to do with CUDA. And the CUDA dgemm should not be used by default for MatrixMatrixMultiply or by `.` or `^`.

I really don't recommend that anyone load the CUDA package anymore. Its supported BLAS functions are mostly slower that the in-memory MKL equivalents on modern machines. If it doesn't serve as a workaround for this example of a dgemm crash then there's no benefit in using it.

The userinfo for LinearAlgebra still shows NAG names, but for all the functions under discussion in this thread the NAG external functions were replaced by Intel MKL equivalents by Maple 18. So I suspect that what's key here is some platform specific behaviour of the MKL dgemm bundled with Maple 2021.

It's useful to know that the multiplication example also crashes for the OP. But it's unfortunate, and surprising, to see a dgemm problem. (Upon my initial reading of this Question I'd expected a problem with the MKL svd function, which has been off-and-on and and occasionally problematic for very large matrices in some lapack & MKL minor releases. And for that routine there is a long-winded workaround to force an older "generic" lapack version of the function. So it's unusual and noteworthy, if instead what we have here is a dgemm matrix-multiplication crash.)

It's usually helpful to try and narrow down the hardware architectecure or operating system specifics on which the crash occurs. Eg, no problem on my i5-7400 running ubuntu 18.04.

I forget whether the following Maple call will work on MS-Windows:
     ssystem("processor full");
But if not then sharing the CPU chipset and OS information might help narrow it down.

Obviously there's no need for workarounds for anyone whose machine doesn't exhibit the problem. For someone who's machine has the problem, I'd be interested in knowing whether the following avoids the dgemm crash. The basic idea is to have the OS environment variable OMP_NUM_THREADS set to 1 in the shell that launches Maple. You'd know it worked if the Maple command,
   getenv(OMP_NUM_THREADS);
returned "1". By default the MKL external libraries intialization functions detect the number of usable cores, when that variable is unset. Setting it to 1 can disable MKL parallelism. For a 2500x2500 Matrix that will make the A^2 example run about 2-3 times slower in wall-clock real-rime on a quad-core machine. But it would be interesting and possibly useful if it prevented the crash.

@MapleUser3214 You didn't state before that you wanted them in an array. That's not quite the same as having two plots printed serially.

Please provide a complete example of the actual kinds of non-empty plots that you want to look better in some array configuration, and the grouping command currently attempted (and which exhibits the problem).

@Kitonum Shorter is the following (and related variants),

    collect(expr, exp, simplify);

But these are techniques appropriate to only a subclass of problems, and are also ad hoc in that the collected terms (here, exp calls) are manually chosen. They don't address the matter of whether desired trig-only simplifications can in general be attained without some measure of factoring/expand that the OP seems to wish to avoid.

@nm You have provided a different example from your original Question.

A similar approach may be take, for `*` as was done for `^`.

restart;

TypeTools:-AddType('specprod',(e,etype,ptype)->type(e,{etype,'`*`'({etype,ptype})}));

type(2*sin(x), 'specprod(specfunc(sin),anything)');

            true

type(sin(x), 'specprod(specfunc(sin),anything)');

            true

This could be used together with the earlier auxiliary type (for `^`) from my Answer. For example,

restart;

# Once these are done, we may add further new types which
# utilize them.
TypeTools:-AddType('specprod',(e,etype,ptype)->type(e,{etype,'`*`'({etype,ptype})}));
TypeTools:-AddType('specpow',(e,etype,ptype)->type(e,{etype,etype^ptype}));

# For example,

TypeTools:-AddType('my_type',
                   'specprod'('specpow'(specfunc({sin,cos}),rational),rational));

type(cos(x), my_type);

            true

type(sin(x)^(3/2), my_type);

            true

type(cos(x)*sin(x), my_type);

            true

If you wanted you could do a similar thing for `+`.

You could call those auxiliary types something else. Perhaps something that didn't look at all like "specfunc" might be less confusing. It's your choice.

@mmcdara I'm not entirely sure I understand where you're going here. Allowing arithmetic operators to be used directly in composing/constructing operators provides for a richer yet also legible syntax. I suppose one could call that a concession of sorts..

@mmcdara I'm afraid that you have missed the main point, and your results are incorrect.

If one is going to work with the result from computing the integral at a symbolic name for upper endpoint then working precision needs to be at least as high as 24 prior to both performing that integration as well as calling plot. Here is one simple way to revise your worksheet, to illustrate my point: PrecisionIssue_ac.mw

However, there are several other interesting aspects, even on top of that corrected approach. There are actually at least three ways to compute the plot with decent accuracy, and the hit-and-miss raising of Digits is only one of them. They include:
1) the fact that exact lowercase symbolic int actually produces accurate results here (using even default-precision float coefficients for the integrand) using exact methods carefully after it detects the floats. That's what happens in my Answer's two suggestions, to use operator form or to delay evaluation (until the passed argument is actually numeric).
2) Forcing purely numeric integration, which happens to work decently accurately here even at its defaults.
3) Raising working precision (eg. Digits>=24) prior to an allowed computation of the integral at a symbolic name for its upper endpoint.

Number 3) is my least favourite, because a better/adequate choice of Digits was found by experiment. (Using shake/evalr/whatever to programmatically ascertain a minimally working precision would take me much longer than trial-and-error...)

Here is my analysis and details which I promised earlier, and examples of the various approaches outlined in the bullet points immediately above.

int_q_details.mw