@Carl Love I believe that the remember table helps performance in the case that plot3d is called with the parameter value varying at inner level.

It removes need to be careful about the plot3d calling sequence, at the expense of memory, etc.

I thought that I'd checked, that both solution approaches provide similar, significant performance improvement. 

@Rouben Rostamian  The key is to have parameter delta change less frequently.

Each time the parameter value is changed then the single dsolve-generated procedure loses knowledge of computation (eg. from initial value up to last t-value requested). Hence it's a disaster to have the parameter be changed on the inner loop, as prior work gets unnecessarily recomputed. Carl's solution is to instantiate and retain/store the procedure for each separate parameter value, while mine is to ensure that the parameter changes in the outer loop.

Another possibility is to construct a Matrix/Array of height values, using independent t and delta, then pass that to plots:-surfdata. That approach can make it more explicitly clear that parameter delta can change in the outer loop, ie. less frequently.

This reminds me that I once started on a rough parallel 3d-surface-plot mechanism, optionally allowing parallelization via Grid or Threads packages. For computationally intensive functions the Grid approach can be effective, which might cover this discussion's example. For computationally lightweight numeric functions using a Compile'd/evalhf'd mechanism that utilizes Threads can be effective. (I think I even prototyped using IterativeMaps and a single iteration, for the latter, since it has some of the machinery already handy. The key is to Compile/evalhf the whole mechanism that populates the Matrix, and not just the function to be plotted.) Maybe I'll find some time to dust all the off and make it more generally workable.

@lcz Yes, I knew that before. That's why I didn't bother suggesting that way.

The mixed dashed/solid effect could be constructed from separate lines. It's very easy to do, but slightly less graceful.

Of course, it's only good for a static rendering, and loses its virtue when rotated -- however that's also true for the textplot placements.

This is why I mentioned the partial dashes in the first place. So far, we don't know whether the OP needs them.

@mclaine I toggle the Mapleprimes editor into Source mode, and then use <pre> ... </pre> tags to get that effect.

Mapleprimes used to support inlining of uploaded worksheets with nice 2D pretty-printing, but unfortunately that's been broken for some months now. (It got fixed for a few days, then broke again.)

@Preben Alsholm The difference in the assignment case is interesting.

restart;
f:=proc() local x; x; end proc:

f() + f() + f() - f();

               2 x

foo := f() + f() + f() - f();

          foo := 3 x - x

restart;
f:=proc() local x; x; end proc:

foo := f() + f() + f() - f();

          foo := 3 x - x

f() + f() + f() - f();

               2 x

There is also some degree of automatic simplification occuring within an additional procedure construction. (There is some match up here, with the behavior given by the OP, which had assignments done at the top level.)

restart;
f := proc() local x; x; end proc:

H := proc() f() + f() + f() - f(); end proc;

    H := proc () 3*f()-f() end proc

H();

            3 x - x

restart;
f := proc() local x; x; end proc:

H := proc() f() - f(); end proc;

    H := proc () f()-f() end proc

H();

              x - x

restart;
f := proc() local x; x; end proc:

H := proc() f() + f() - f() + f(); end proc;

      H := proc () 2*f() end proc

H();

               2 x

@mmcdara Thanks, but I had already run in the debugger, and already knew exactly how many times f() gets evaluated in these cases. The debugger cannot help much more than that, here.

An issue with these mixed cases (positive and negative coefficients) is inconsistency between cases and between what happens after wrapping in uneval quotes (relatively speaking). Some amount of automatic simplification occurs. 

The debugger can only show what happens when the evaluation occurs, and that happens only *after* any automatic simplification occurs.

@Preben Alsholm As far as I know (so far) that too is a result of some automatic simplification.

See also the cited example,  f()+f()+f()-f()

The behavior is more complicated in these mixed cases. [edited] Earlier I wrote, "...affected by the order of terms in the sum" but I am not sure. See later comments below.

@mmcdara In general, symbolic integration for expressions containing floating-point coefficients is a not a great idea. This general problem relates to broader areas of symbolic computation, and not just integration in statistics or DEs.

You have omitted to show how M and L figure in your conditions, following "such that".

The wrong use of the `assume` command and the handling of names inside procedure S is the main cause of the problem.

Using the deprecated linalg package, setting Digits to 2, and use of `%` inside procedure S are also poor choices.

@Kitonum The conjugate is there. I suspect that the later assumptions were one way for the OP to deal with it.

If that subterm were not conjugated then it could be done more directly (ie. without substitutions or supplying extensions to factor) as,

   simplify(evala(Expr), symbolic)

or,

   simplify(radnormal(Expr), symbolic)

or even,

   simplify(expand(rationalize(Expr)), symbolic)

@Kitonum Did you omit the conjugate term in entering your version of the expression?

@lime I know you are trying to optimize performance time. That's why I suggested reducing the accuracy tolerance. But no accuracy value optimizes all possible examples, so fishing for a value that works for some unknown integrand is difficult.

Also, it may well be that the computation of the integrand could be sped up. But that's really hard to do without knowing what it is. Is it evalhf'able?

Specifying some methods can save time by avoiding unnecessary checks for discontinuties. (Another way is using so-called operator form.)

This is a complicated and general query, without the integrand being known.

[edit] I expect that you've already tried exact integration.

There is often very little difference between properly using evalf(Sum(...)) and add(evalf(...) for a small finite number of terms. I doubt that by itself would help in a significant way.

I'd suggest instead to focus first on the performance of individual calls to numerically compute the integrand -- especially under evalhf, but also more generally. And, as originally mentioned, loosening the tolerance,  and forcing a method or using operator form.

It usually helps to explicitly supply the integrand. You can upload and attach a link to a worksheet here by using the green up-arrow in the Mapleprimes editor.

Some performance improvement might be attained by passing options such as, say,
    epsilon=1e-4, method=_d01ajc
to the Int call.

This Answer makes claims about the original results from the OP. But actually it is this Answer which is wrong.

I don't understand why some people upvoted this Answer.

There is no justification provided for throwing away additional digits of precision in the data. If one keeps the full digits provided and accomodates the numeric error (eg. roundoff, loss or precision) then by raising the working precision the original plot is explained.

Conversion to rational without the exact option will throw away information here, without justification. Conversion to rationals is not actually necessary here, if Digits is raised adequately. But too low a wokring precision for such rational conversion gets highly inaccurate and wrong results for this example.

prec_examp.mw