You might look at some of the approaches in these responses to these old postings: 1, 2, 3

note: Considerable portions of those approaches can now be done more straightforwardly, ie. in modern Maple. For example, using the background option of the plot command, using convert(plot(...),Image) rather than exporting a plot to an external temp file, etc. When I first used an image here, in the year 2012, even the ImageTools:-Embed command did not yet exist, so I had to first insert a Label Component manually.

How are you currently exporting your plots, for use in LaTeX?

What specific steps (mouse+GUI, or programmatic commands) are you using? What target file format? What resolution?

@paulmcquad For your examples, you can also use mmcdara's straightforward Answer.

@Kitonum I gave it a vote, and marked it as a favorite.

@erik10 Yes, the integrand has a singularity where z^2/(p^2*v__lin(z)^2) - 1, and is complex-valued for z making that expression negative.

For your positive data, at z=p*v__lin(z) .

Seismology_acc_ev_ac_ev_ac.mw

@erik10 So what is happening in your new sheet is that you are passing x(r) and y(r) to the plot command, and those evaluate to expressions involving theta(r).

But your procedure assigned to theta is not prepared to be called with symbolic r (ie. not an actual number). If you were to issue just theta(r) as a separate statement then you'd see the same error. So, those function calls passed to plot are evaluating prematurely.

One way around that is to use the so-called operator-form calling sequnce of the plot command, eg,
    plot([x, y, r__0 .. R], ...options)

Another way is to instead construct the theta procedure so that it returns unevaluated when called with a symbolic argument such as just a name. You could do that using MakeProcedure (unapply) with its numeric option, as I showed before for your Z__spline, etc. Or you could write it all out, which I do in the attachment.

I also add the numeric integration option epsilon=1e-5 here, which relaxes the accuracy tolerance for the theta computation (the outer integration). That makes your nested numeric integration faster, as I mentioned earlier that I suspected. It's now fast enough that you can get smooth adaptive plotting without too long a wait.

Seismology_acc_ev_ac.mw

@erik10 The adaptive=false option instructs the plot command to not do adaptive sampling of the domain. It makes it use an evenly spaced number of values in the domain.

The numpoints=N option instructs the plot command to use N values in the domain, initially. If adaptive=false then that means N evenly spread values, and otherwise it means N values initially in whichever adaptive sampling scheme it's using.

In other words, if one wants only N data points, evenly spaced, then one can supply,
   adaptive=false, numpoints=N
This is a way to know exactly how many evaluations of the function will be done, so is handy when it might be time consuming. When going this route I experiment to find the smallest number of values which make the curve appear smooth "enough". (Sometimes I start with N=3, as sanity-check experiment that things are working at all...)

Of course, an alternative is to use seq or Vector to get a specific number of evenly spaced data points. But I prefer to use the above options since if everything goes well or if I speed it up then I can just remove those options.

@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.