@planetmknzm Given two points on the unit sphere, the central angle (angle at the origin) between them is the arccos of their dot product. If the points are represented as lists P and Q, that's 

arccos(<P>^%T  . <Q>)

There is no picture attached to your Question.

If you have a Maple worksheet, it'd be more helpful if you uploaded that. Use the green up arrow on the toolbar.

@planetmknzm Thank you. I finished the explanation. Let me know if you have any questions.

@planetmknzm 

The geometric simplicity of the surfaces being plotted allows some of the commands that I used to become deceptively simple. The first such simplicity is:

• Given a point P on a sphere centered at the origin, if P is viewed as a vector v, then v is a normal vector to the sphere at P and hence also a normal vector to the plane tangent to the sphere at P.

I'll detail the other simplicity later.

Let's look at the details of an animate command:

plots:-animate(PP, [A(t)], t= a..b, frames= n)

• PP is a procedure that returns a plot.
• A(t) is the arguments to PP parameterized by t.
• a and b are real numbers, the range of values for t.
• frames is a keyword, i.e., it appears literally.
• n is a positive integer, the number of evenly spaced values of t that'll be used.

My animate command was 

plots:-animate(Plane, [C(t)], t= -Pi..Pi, frames= 100, background= Static)

C(t) is my command that generates points on the curve on the sphere. So, 100 such points are passed to Plane. Each point is a list of three real numbers.

Now let's look at Plane:

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

A list can be turned into a column vector by enclosing it in angle brackets: <P>. The transpose operator ^%T converts column vectors into row vectors and vice versa. So, v is a row vector, the normal vector of the plane. The distinction between lists, column vectors, and row vectors is of very little mathematical significance; these conversions are just for syntactic reasons.

When NullSpace is passed a row vector v, it returns a basis for the subspace of vectors that are orthogonal to v. In this case, that means the vectors in the tangent plane. The basis is returned as a set of column vectors. I want them to be row vectors for the next command. By adding the elementwise symbol ~ to the transpose operator ^%T (so that the operator is now ^~%T), the entire set is converted to row vectors. So B is now a set of 2 row vectors orthogonal to v.

Now here is the second hidden simplifying detail:

• The NullSpace command returns an orthonormal basis, that is, the basis vectors are orthogonal to each other and they have length 1.

So, adding and subtracting the basis vectors to the point of tangency gives the four corners of a square. That square has diagonal length 2, thus side length sqrt(2). This seemed to be a good size for the animated squares, but you could change it by multiplying the vectors in B by a scaling factor.

The command polygonplot3d accepts a three-column matrix each row of which represents a vertex of the polygon. (The rows are viewed as points, not vectors.) The angle bracket syntax <..., ..., ..., ...> stacks the row vectors into a column, which is the same thing as a matrix.

@Ronan In the following example, note how the invisible function application operator distributes over the comma sequencing operator between the f and g:

(f,g)(x,y);
                        f(x, y), g(x, y)

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