Your attached worksheet Q20220425.mw is empty.

Great. After I mentioned parallelization of Jac1 I realized that it was generating a Matrix with storage=sparse and that likely it'd be more key to construct it via compressed form (a triple of vectors). Your second worksheet appears to be working along those lines, so that's a place to focus. For your earlier sheet I'd noticed that the cost of construction by Jac1 was even more important than that of the lienar-solving as NN grew.

For your first worksheet my machine takes 25.7sec for NN=100 (dim of j00=10404, with WEN03).
After removing overhead of calling UMFPACK's direct solver, my machine takes 22.7sec for the same NN=100. Phasefield2DBaseCodeParametric_ac0.mw
Running your 2nd worksheet my machine takes 11.3sec for a similar ntot=100.

I would like to look closely at construction of the sparse Matrix ("jacobian"), especially via direct compressed formats.

I'd like to look at replacing both hw_PardisoFactor and hw_SpUMFPACK_MatFactor with some custom wrapper revisions that might allow reusing preconditioner/"factorizer" construction.

I have a hazy recollection of once reading that a sparse iterated approach might begin to beat a sparse direct approach (say, for this kind of block-banded solve that comes up so often in these kinds of pde problems) as the matrix dimension grows above about 10k-by-10k. It would be interesting to investigate any such possible cross-over. But that would of course require using the very same "jacobian" construction for both. And ideally that'd also need optimal "re-use" for both, as mentioned.

I might get some time, starting tomorrow evening.

note: Out of curiosity, your 2nd worksheet was apparently saved last by you in Maple 2021.1, and you use Linux IIRC, and it contains a define_external to hw_PardisoFactor. Did you run it in 2022.0 or actually in 2021.1?

note: If I have to compile any new custom wrapper would it be ok (say, to start with) if I supplied a linked .so file (Linux) rather than C source?

The possibility for improvement looks likely.

In the best case scenario the common pattern of sparse Matrix j00 might be leveraged, by reutilizing a "symbolic" preconditioner generated by UMFPACK. But this may require effort, as Maple's pre-built external wrapper function for that initial stage returns only handles (to objects in memory) which might not be thus re-usable. A new hand-crafted wrapper might be constructed, utilizing umfpack_di_symbolic once and then umfpack_di_numeric for each new Matrix.

That could be quite beneficial if the three linear-solving computations altogether consume about 2/3 of the time for each iteration at that NN=100.

On the other hand, a certain amount of Maple overhead might be cut out, though the savings for that might not amount to so much --> perhaps 20% of the total time, at that size.

But if the problem size is increased then it looks as if some of the other (Compiled, custom) procedures such as Jac1 might benefit relatively more from improvement -- possibly including parallelization.

So here are three ideas: re-usable factorization, reduction of Library call overhead, and parallelization of the custom procedures. I can see definite paths for the latter two.

Please add any additional details (or close followup queries) here, instead of spawning separate new Question threads on this.

What have you managed to accomplish so far, for this homework?

Please don't spawn a separate Question thread for the same kind of query. Instead, please put your additional details or close followup example(s) here, in Comments/Replies.

Are you sure about your (plaintext) code, which you gave as, eg.:

find_vars_in_proc(f :: procedure, $)
  return {op(2, eval(f))};
end proc;

Also, there doesn't seem to be much merit in creating a procedure whose actions could just as easily be accomplished directly (without recourse to a procedure).

There are several syntax mistakes in your formulas. There significant details missing in your Question.

Most of the syntax mistake appear to be with respect to multiplication. For example, does gamma(nu-1) without any space of explicit multiplication sign mean gamma*(nu-1) ?

You should fix all those syntax mistakes first.

In Maple the lowercase gamma is a predefined constant, ie. evalf(gamma)=.5772156649. If instead you intend it as some other (free) parameter then you could use a different name instead.

You have not explained which names are the main plotting variables and which are the free (explorable) parameters.

You have not explained what kind of plot you want. Do you want a pair of 3-D surface plots for each of E[1] and E[2], as functions of say mu and nu? Or do you want a 2-D parametric plot of E[1] versus E[2] as a function of rho? Or do you want something else? There are several other reasonable possibilities, and it's not possible to know for certain which is your goal, unless you explain it clearly.

You should explain your intent here, clearly and explicitly.

@Carl Love It was fixed between Maple 2015.2 and Maple 2016.2.

@Arastas The trailing tilde makes the command map across entries of a Matrix, or list, etc, in an elementwise manner.

Some commands automatically get applied across elements of such structures, but most do not.

@etian2 It works for me. In 2D Input mode I type in,

  a := 5

but before I press <Enter> I select the "a", make it bold, then right-click to make it an atomic identifier. Then I press <Enter>.

Then in a new statement I type "a", repeat the process to make is bold and an atomic identifier. And after I press <Enter> the evaluated output shows as 5.

Download bold_assignment.mw

Why do you think that your procedure is thread-safe, given that it calls signum?

I don't see signum listed on the Help page for Topic threadsafe.

I don't know whether that is part cause of your problem.

It does not make sense to say that you cannot upload a Maple 11 worksheet because "Too complex to insert the Maple 11 version because I have  a complicated function  and further on".

The integral you claimed to work in Maple 11 is short. It either works, by itself and simple, or it doesn't.

@Preben Alsholm I get the following:

restart;

kernelopts(version);

    Maple 11.02, X86 64 LINUX, Nov 9 2007 Build ID 330022

with(IntegrationTools):
with(DEtools):

fx:=exp(-sqrt(x^2 + c)):

integral:=int(fx,x) assuming positive;

      int(exp(-(x^2+c)^(1/2)),x)

value(integral);

      int(exp(-(x^2+c)^(1/2)),x)

int(fx,x=xl..xu,method = _DEFAULT) assuming (x^2+c)>0,x>0,c>0;

              2     1/2
  int(exp(-(x  + c)   ), x = xl .. xu, method = _DEFAULT)

int(fx,x=xl..xu) assuming (x^2+c)>0,x>0,c>0;

   int(exp(-(x^2+c)^(1/2)),x = xl .. xu)

It is poor etiquette and unhelpful to show only an image, instead of actual code to reproduce.

You could lprint the expression, and show us that as plaintext, or you could upload and attach a worksheet using the green up-arrow in the Mapleprimes editor.