@erik10 You could use the MakeFunction (unapply) approach that I'd shown above:

Seismology_ac.mw

Or you could use an operator-form for your new integrand, as variantion on what I'd shown before for that kind of idea. (This variant now gains a slight inefficiency of a separate procedure construction for each call to theta, though in practice that can be relatively reasonable if the numeric quadrature is itself expensive enough.)

Seismology_acc.mw

[edit] I suspect that the final plot (of theta) might be computed faster if the epsilon tolerance were loosened on one of the (nested) numeric integrations. But for that you'd need to give up the pretty 2D integral look. I left it with non-adaptive plotting and a prettier integral, in these two worksheet attachments, and not too much of a wait.

ps. Often it's more useful to supply the background example up front.

pps. It's a good question, though, and has come up before. The "evalf-Int" routines fail to deal with the constructed interpolating function (an object) in as versatile a way as they do regular appliable procedures. I should make a(nother) bug report on it.

@salim-barzani If you have similar problems often then you might consider using 1D Input instead (see red plaintext code in my worksheet).

You can set that as default as a GUI Preference, or you could toggle between 1D/2D modes on the fly with the F5 key.

@segfault We gave different ways that produce the same part for your asympt example, ie. the O term.

I wrapped my first way with order and I wrapped my second way with degree@op, as programmatic followup since you mentioned both the O term as well as the order itself.

You can compare:

s:=asympt(x/(1-x-x^2),x);

s-convert(s,polynom);

series(s,x);

indets(s,specfunc(O))[1];

# or,
indets(s,specfunc(O))[];

@Ronan If we did just Metric:=M then some change to the original Matrix would also be seen when subsequently accessing the stored local Metric.

But such a change to the original Matrix would not be accompanied by a recomputation of the inverse. So the stored inverse could then be wrong.

You wouldn't want that behaviour, as in this modified worksheet, I suspect:
2024-12-30_Q_Module_Global_and_Local_oops.mw

@sursumCorda You can try wrapping that part in an eval call, within the proc to be compiled. Eg.,
   eval(`if`(isprime(j),1,0))
The Compile command may issue a Warning about isprime not being recognized, but in this example it seems to work ok.

ps. The evalhf mechanism also allows for this kind of escape (back to Maple proper).

@Axel Vogt That's a fun idea.

That inner loop could be exited as soon as f ever became zero. Eg,

for k to j-1 do
            f := f*k*k;
            f := f mod j;
            if f=0 then k:=j-1; next; end if;
        end do;

That brings the cJ(40) timing down from 17.4sec to 4.7sec on my machine.

axel_J_ac.mw

Naturally, this is still trying to adhere to the OP's formula using mod, not some faster equivalent.

@nm You wrote that you are making relative paths, but that's not true. Your given cat calls make absolute paths.

And the extra step of applying cat is not necessary, since you had already set currentdir. That concatenation entirely misses the point and purpose of using currentdir.

(If absolute paths were wanted you could also just assign the prefix to any name.)

@vv You may well already know this, but pardon me if I mention that one can extract a value using sol(last) in your code.

Eg. colibri-vv_ac.mw

ps. I'm reminded of this old post, where I was just having fun with dsolve's events option.

@Rouben Rostamian  When n=4 and j=16 then the compiled code tries to compute (15!)^2 which is too big for a 64bit integer in compiled C, and it throws that error.

So one might ask: how to compute ((j-1)!)^2 mod j without having to actually compute ((j-1)!)^2 explicitly.  Handling only the square (via inert `mod` stuff`) is not enough to go much further; the factorial is a sticky wicket.

@ecterrab If I run the same worksheet as I did before, using Maple 2024.1, then some of those diff calls take a while, and then throw an error.

So there was certainly some point release at which it did not work ok in Maple 2024.

I mention all this so that the OP of this Question understands that loading the Physics package can indeed make a difference here.

I leave it to the OP to decide whether to go with Maple 2024.2, or not load Physics, or install the latest add-on Physics update package, etc.

Download why_ac_M2024.1_oops.mw

@salim-barzani It works ok in my Maple 2024.1 without Physics loaded.

Download why_ac_M2024.1.mw

@nm I suggest using Grid instead of Threads for dsolve computations that you want to run in parallel.

What happens if you omit the with(Physics) line?

The limit and int commands are not Thread-safe, in general, and neither is dsolve.

May I ask, did you file a Software Change Request (ie. bug report) against this, either by form, direct message to Tech Support, or via the Beta forum?