@ivar_kjelberg You wrote:

• I got the impression that "init" was the old way, and "ModuleLoad" the more recent implementation.

It is true that init is older than ModuleLoad, but, as I said, they are independent. They do different things, and ModuleLoad is not a replacement for init; nor has init been labeled as deprecated. They are different because ModuleLoad() is triggered by the module (which is not necessarily a package) being extracted from an archive, whereas init() is specific to package modules and is triggered specifically by the with command.

You did ask specifically about the with command, so, to directly Answer your title Question:

• Does a ModuleLoad proc() run with the "with()" command?

Direct Answer: No, ModuleLoad() is not necessarily invoked by the with command, but init() is.

• I even cannot find any documentation mentionning the init for Maple 2022.

Yes, it's hard to find. It's the 11^th bullet point in the Description section of help page ?with:

• Prior to arranging for the names exported by a package to be available at the top level, with executes either or both of two procedures that have, for historical reasons, come to be known as the system level initialization and user level initialization routines for a package, provided that they exist. A system level initialization routine is any procedure that is assigned to the name P/init, where P is the name of the package. It is called with no arguments. The user level initialization routine is a package export called init. It is called (if it exists) with no arguments after first calling the system level initialization routine. The system level initialization procedure is not called if the name P/Initialized has the value true, and this name is assigned the value true once the routine has been successfully executed. This allows you to put any once-per-session setup that needs to be done in the system level initialization, and anything that must be run each time with is called on your package can be performed in the user level initialization.

@Nikol Clearly there is a bug here. I'm not recommending that you use expand; that's too ad hoc to be generally useful. However, you're rightfully curious about why expand works. It's because it converts exp(-T) to 1/exp(T) and exp(-2*T) to 1/exp(T)^2. The expansion of the polynomial parts is irrelevant. Letting Y= exp(T) and clearing denominators allows the equation to be rewritten as a polynomial that is quadratic in two variables:

AY:= (20 + 20*T + 2*T*(T+1))*Y - 10 - (2*T + 10.0)*Y^2 = 0:

solve({AY, Y=exp(T)});
             {T = 1.411454823, Y = 4.101918631}, {T = 0., Y = 1.}

Once again, this is too ad hoc, so I'm not recommending this in general. I'm simply presenting it as a curiosity.

@rlopez The connect-the-dots phenomenon that you mentioned is often an issue with the plot command. (In Maple 2022, the new option setting adaptive= geometric, which is now the default, has largely corrected this.) But the case at hand uses implicitplot, and inspection of the plot structures reveals that that phenomenon is not the issue here. Two of the shown plots have a vertical line. (They are truly vertical lines in these cases, unlike the nearly vertical lines produced by connect-the-dots.) Looking at the plot structures for the x= -6..9 and x= -6..11 cases shows that they each contain three explicitly separate CURVES, one of which is (in each case) a vertical line composed of 201 points. In the x= -6..9 case, that vertical line is x = 1.96569296691827; in the x= -6..11 case, it's x= 2.05767071950135.

I don't know why these vertical lines are drawn. Their points are not even close to being zeros of the input expression y - x/(2-x). But it's clear that it's the high-level Maple-language plotting command that's making the error, not the plot renderer.

@fkohlhepp Option numeric is for int, not Int. Your original Question didn't show your integration commands, so I didn't know that you were already using evalf(Int(...)), which is pretty much equivalent to int(..., numeric) anyway. So, you're already using numeric integration, so switching to int(..., numeric) should "work" superficially in the sense of not giving an error, but it won't help with the computation time. So, there's no point in you using numeric.

Given that you're already using numeric integration, the slowness may be due to one of these causes:

1. Singularities: The integrand has singularities within the interval of integration (essentially, values of r where division-by-zero happens or where the integrand is not differentiable). These might still be integrable, but the algorithms to do it numerically are more complicated.
2. Loss of precision: The integrand is complicated enough that the integrator is having trouble evaluating it numerically at sufficiently high precision to get the required precision for the final integral.
3. Massively complicated: The integrand is simply so complicated that each numeric evaluation takes a long time not due to precision but simply due to the tens-of-thousands of steps required.
4. Inadvertent symbolic computations: Due to inadvertent coding of your procedure trq__s, it may be needlessly redoing symbolic computations for each numeric call. One of the most-common cases of this is a procedure which computes a derivative before evaluating it numerically.

Do you have an idea that one of those is the culprit? Showing an example plot of the integrand over the interval of integration and recording the time required to generate that plot will help.

Using some or all of these might help:

1. Put the Re(...) inside the integral.
2. Add the option epsilon= 0.5e-5 to the integration command. This will reduce the precision (precisely, increase the maximum allowed relative error) (of the final result, not of the internal computations) to 5 digits. That number can be adjusted. Even if this level of precision is unacceptable for your purpose, doing this just as a test will help to track down the source of the slowness. 
3. Add one (and only one) of these options to the integration command: method= _d01ajc or method= _d01akc. These are high-speed externally compiled numeric integrators from the Numerical Algorithms Group (NAG). 

@C_R The ' 'arctan' '(algebraic) is a structured type which is essentially used as the 2nd argument of an indets command. That indets command searches for function calls with function name arctan having exactly one argument of type algebraic (see ?type,algebraic).

Since arctan is an actual executable procedure (not just a function), specifying arctan(algebraic) would cause arctan to execute, which is potentially undesirable. (It actually wouldn't be a problem in this particular case because it would return literally arctan(algebraic) unevaluated. But why waste the execution time required for that?). The quotes prevent that execution. One pair lets the expression pass through rcurry, and the second lets it through evalindets.

To understand what the builtin evalindets does, it's useful to rewrite it in top-level Maple. When it has three arguments and no options, it's essentially the same as this:

EvalIndets:= proc(e, T::type, f)
local r:= e, k, Ind:= sort([indets(e, T)[]], 'key'= length);
    for k to nops(Ind) do
        (r,Ind):= eval([r, Ind], Ind[k]= f(Ind[k]))[]
    od;
    r
end proc
:

Now we can trace its execution for educational purposes:

`convert/atan2`:= rcurry(
    EvalIndets, ''arctan''(algebraic), arctan@op@[numer,denom]@op
):
eq:= phi = arctan(A/B) - arctan(arctan(A/B)/B):
trace(EvalIndets):
convert(eq, atan2);

You should study the output of that. It's only 11 lines.

@Stretto I don't know if that's possible via Maple commands. Essentially, you would like to click somewhere and have Maple respond with the 2D-coordinates of where you clicked. (And that's a very natural thing to want; I've thought about it for years.) If there were a formal mechanism for this, I think that it'd be in the DocumentTools:-Actions subpackage, but it isn't.

The functionality must exist within the GUI at some level, because you can do this: Using the context menu for a 2D-plot, select Probe Info, then Copy Point Data (to clipboard). I just did that (on some random plot), and I got this:

-4.421
2.172

In other words, I simply pasted my clipboard to the MaplePrimes editor rather than literally typing those numbers. So I suspect that there's some hack possible to get those numbers programmatically. Acer @acer might know how. If anyone here knows, it'd be him.

As acer just hinted at, arctan(y/x) and arctan(y,x) are not necessarily equal. For x<0 and real y, these expressions will differ by Pi; however, for x>0 and real y, they'll be equal. I was only providing a simpler means to do what you were already doing by more-complicated methods in your worksheet. The information about angles in the 2nd and 3rd quadrants has already been irretrievably lost by the use of the one-argument form.

@brownr Sorry, I misread the Question. Yes, you and the OP are correct that Maple Flow is not listed under that Products tab. And, of course, it should be listed there.

@nm Several numerical calculations that I've done (with the ODE with initial condition y(0)=2) suggest that using the process that I described produces a correct solution to your original ODE and that the other solution (the one with hypergeom) is complete garbage that doesn't even come close to satisfying its own initial condiiton (y(0)=2)! 

For example:

restart:
ode1:= diff(y(x),x) = x*(y(x)^2-1)^(2/3):
ode2:= diff(y(x),x)^3 = x^3*(y(x)^2-1)^2:
ic:= y(0)=2:
sol1:= dsolve({ode1, ic}); #hypergeom solution
evalf(eval(rhs(sol1), x=0)); #verify initial condition
               -1.005937890 - 4.630275440*10^(-14)*I

#Expected 2., possibly with some rounding error.

sol2:= dsolve({ode2, ic}); #WeierstrassPPrime solution
#Check in original ODE:
res:= odetest(sol2, {ode1, ic});
#That residual is extremely complicated, but plotting suggests that it's
#identically 0, except for some "noise" that's totally expected due to
#rounding error when plotting such a long and messy expression. 
plot(res, x= -2..2);

So, to directly Answer your initial Question, to wit:

• Any one knows of a way to have Maple verify that this solution of the ODE is correct?

Answer: It can't be verified as correct because it's actually not correct.

This is a blind guess; I have no way of testing this: I suspect that the 0 is the operating system error code for a successful (i.e., without error) run of an external program. If you were using Maple instead of MapleSim, an analogous command would be ssystem, whose return value is a 2-element list, the 1st element being that error code, and the 2nd element being the actual data returned by the external program.

If this analogy isn't just coincidental, then you need to somehow extract the 2nd element of the returned data.

@lcz And here is another use of RelationGraph: a procedure for GraphPower with 3 methods:

1. an algorithm that I just invented called unions;
2. a method based on the AllPairsDistance matrix and RelationGraph;
3. Maple's own GraphTheory:-GraphPower.

My unions algorithm in several times faster than the others for small values of k.

RelationGraph:= (f, V::list({integer, symbol, string}))->
local n:= {$nops(V)}, i, F:= i-> j-> f(V[i], V[j], _rest);
    GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
:

GraphPower:= proc(
    G::Graph, 
    k::posint,
    {method::identical("default", "distance", "unions"):= ""},
    $
)
local V:= GraphTheory:-Vertices(G), n:= [$nops(V)], N, L;
    if method = "unions" then
        L:= rtable(op~((N:= op(4,G))));
        to k-1 do N:= N union~ (n-> {seq}(L[[n[]]]))~(N) od;
        GraphTheory:-Graph(V, N minus~ `{}`~(n))
    elif method = "distance" then
        N:= GraphTheory:-AllPairsDistance(G);
        L:= table(V=~ n);
        RelationGraph((u,v)-> u<>v and N[L[u],L[v]] <= k, V)
    else
        GraphTheory:-GraphPower(G,k)
    fi    
end proc
:

And here is a usage test comparing all methods:

GT:= GraphTheory:  RG:= GT:-RandomGraphs:
G:= RG:-RandomGraph(2^10, 2^11):
gc(): CodeTools:-Usage(GraphPower(G, 5, method= "default")):
memory used=71.82MiB, alloc change=6.01MiB, 
cpu time=1.16s, real time=1.15s, gc time=0ns

E1:= GT:-Edges(%):  nops(%);
                             315735

gc(): CodeTools:-Usage(GraphPower(G, 5, method= "unions")):
memory used=60.46MiB, alloc change=0 bytes,
cpu time=281.00ms, real time=280.00ms, gc time=0ns

E2:= GT:-Edges(%):  evalb(E1=E2);
                              true

gc(): CodeTools:-Usage(GraphPower(G, 5, method= "distance")):
memory used=240.44MiB, alloc change=6.01MiB, 
cpu time=3.69s, real time=2.99s, gc time=812.50ms

E3:= GT:-Edges(%):  evalb(E1=E3);
                              true

@lcz You wrote:

• I am probably more interested in the panning implementation of RelationGraph in maple. Maybe it's not always the most effective. But it can reflect the degree of abstraction of certain graph operations. And not just some graph operations in isolation.

I agree on all points. A procedure for RelationGraph is quite easy; here it is:

RelationGraph:= (f, V::list({integer, symbol, string}))->
local n:= {$nops(V)}, i, F:= i-> j-> f(V[i],V[j]);
    GraphTheory:-Graph(V, [seq](select(F(i), n, _rest), i= n))
:

And here is a rewrite of LexProd which offers a choice of two methods: one using RelationGraph and the other the same as the algorithm in my previous Reply:

LexProd:= proc(
    G::seq(Graph),
    {
        method::identical("default", "relation"):= "default",
        vertex_format::{string, symbol}:= "(%a,%a)"
    },
    $
)
local
    LexProd2:= proc(G::Graph, H::Graph, $)
    local 
        (Vg,Vh):= op(GT:-Vertices~([G,H])),
        (ng,nh):= op(`..`~(1, nops~([Vg,Vh]))),
        (Ng,Nh):= op(op~(4, [G,H])),
        i, j, J, k, K,
        P:= [seq](seq([i,j], j= nh), i= ng),
        Vgh:= (curry(nprintf, vertex_format)@((i,j)-> (Vg[i],Vh[j]))@op)~(P)
    ;
        if method="relation" then
            K:= subs(_K= op~(table(Vgh=~ P)), v-> _K[v]);
            RelationGraph(
                proc(a, b, $)
                local (u,v):= K(a), (x,y):= K(b); 
                    x in Ng[u] or u=x and y in Nh[v]
                end proc,
                Vgh
            )
        else
            k:= 0;  K:= table((p-> op(p)= ++k)~(P));
            GraphTheory:-Graph(
                Vgh,
                [seq](
                    {seq(seq(K[k,J], k= Ng[j[1]]), J= nh), seq(K[j[1],k], k= Nh[j[2]])},
                    j= P             
                )  
            )
        fi
    end proc,
    n:= nops([G]);
;
    if n=0 then error "no graphs given" elif n=1 then G[1] else foldl(LexProd2, G) fi
end proc
:

And here's an example verifying that the two methods give identical results:

GT:= GraphTheory:
(G,H):= GT:-Graph~([`$`]~([3,4]), [{{1,2}}, {{1,3},{2,4}}])[];
 G, H := Graph 38: an undirected unweighted graph with 3 vertices and 1 edge(s),
     Graph 39: an undirected unweighted graph with 4 vertices and 2 edge(s)

GH1:= LexProd(G,H);
    GH1 := Graph 40: an undirected unweighted graph with 12 vertices and 22 edge(s)

GH2:= LexProd(G, H, method= "relation");
    GH2 := Graph 42: an undirected unweighted graph with 12 vertices and 22 edge(s)

evalb(GT:-Edges(GH1) = GT:-Edges(GH2));
                              true

@Earl You wrote:

• I've spent some time studying your latest replay. It's difficult for me to understand so I'll return to it later.

I'd be happy to answer any questions about it, especially if it's not a large number of questions asked all at once. (However, a large number of questions spread out over several sessions is not a problem.)

@Joe Riel @AmirHosein Sadeghimanesh

It should be pointed out, for the OP's benefit, that

1. Joe's workaround works because MutableSet has a procedure named ModuleIterator;
2. the original command doesn't work because MutableSet doesn't have a procedure named ModuleType;
3. the workaround in my Answer below works because MutableSet does have a procedure named convert which accepts the keyword set as its 2nd argument.

In summary, what works is determined entirely by an object's defining module having specifically named procedures (or "methods" if you want to call procedures that).

@planetmknzm

Yesterday I saw that you asked (essentially) how to compose an animation consisting of some parts that move and some that don't. That Question got deleted, I don't know if by you or some Moderator. Anyway, I'll Answer that Question here.

Here is a program that plots a curve on the unit sphere and animates a plane tangent to the sphere whose point of tangency moves along the curve. As you can see, this can be done much more simply than your attempts. My computer uses 0.17 seconds to produce this 100-frame animation.

restart:
Sph:= (theta,phi)-> [(cos,sin)(theta)*~sin(phi), cos(phi)]:
C:= Sph@op@[theta, phi]:

#Example curve:
theta:= t-> t+Pi:  phi:= t-> t^2/Pi:

#Parts of the animation that don't move:
Static:= plots:-display(
    plottools:-sphere(), 
    plots:-spacecurve(C(t), t= -Pi..Pi, thickness= 4, color= blue)
):
#Plot a plane tangent to sphere at a given point P: 
Plane:= proc(P)
local v:= <P>^%T, B:= LinearAlgebra:-NullSpace(v)^~%T;
    plots:-polygonplot3d(<v+B[1], v+B[2], v-B[1], v-B[2]>)
end proc
: 
plots:-animate(Plane, [C(t)], t= -Pi..Pi, frames= 100, background= Static);