@Nicole Sharp The end user need not know about the local file, in the sense that the file can be downloaded and immediately read every time. (The intermediate file stage could be avoided if the read statement had a source=direct option like import does, which means treat the string as file contents.)

You can take the string, use StringTools:-Split to split it at each endline "\n" character and them parse each line. Parse only does one expression or statement, and then you have to handle the statements differently from the expressions. In the end writing an intermediate file seemed to me the easier option.

(The Import statement is handling the .mpl file in just the same way as if it is reading a local file. It defaults to output=string.)

@WA573 Use

interface(showassumed=0):

to remove the ~

@C_R Yes, equatorial poles are harder with spherical coordinates. If you use phi>Pi/2 and have four regions of theta then you can have each octant colored separately; the equatorial and north/south poles are then at the points where the colors intersect. You could add something like

pointplot3d([[0,0,1],[0,0,-1],[0,1,0],[0,-1,0],[1,0,0],[-1,0,0]],
  symbol=solidsphere,symbolsize=10,color=black):

though there are little bumps then.

An equatorial band is possible with a small phi range and north-south bands with a small theta range, though the north-south ones narrow at the poles.

@Anthrazit  Thanks for the complete information. To make a nested table, you need

WhateverYouNeed["calculations"] := eval(calculations);

Edit: my earlier explanation is too simplistic. In one step:

WhateverYouNeed:=table(["calculations"=
                   table(["structure"=
                     table(["connection"=
                       table(["Cutright1" = "false", "graindirection1" = 90*Unit(arcdeg)
                       ])
                     ])
                   ])
                 ]);

TableDefinition.mw

@WA573 You needed nested seq, and I changed subs to eval, but the calculated values are complex and so you get an empty plot.

plot_dis.mw

@Zeineb I don't see any Maple errors, just a warning that i was not declared in A, which I fixed. You haven't written the integral as an integral with respect to psi, not t - just cancel "dt: in (dpsi/ dt)*dt = dpsi as I did above. I'm not sure what you are trying to do but you still have A's as a function of y's and y's as a function of A, which is confusing.

last.mw

@nm I agree. Indeed, odeadvisor has the option of testing other types, which you would not need if if was intending to return all types.

DEtools:-odeadvisor(ode,y(x),[homogeneous]);

gives

@WA573 In (6)-(8) I only see bars over expressions that contain eta, and are therefore potentially non-real; I don't see any bars directly over n or m alone.

@Zeineb Changing int to Int in Frac_D shows why the integrals are being evaluated as zero. Notice that A[0](x) works as expected but A[0](t) does not. This can be fixed by local t (in red), though this may be a side issue. But the code is confusing because the y[i] are used in the definition of A[n] but are defined later, and I lost track of locals and globals. I suggest you create the A[n] and y[i] one by one  in the right order, and only when it works put it into one loop.

compute_integral.mw

@Zeineb Yes, Maple can do Axel's integral given those assumptions. Simplify (after dividing by GAMMA(alpha)) then gives the form you want.

gives (in Maple 2023):

@nm I agree a better simplification would help here, and I think that was @acer's point - it didn't help. In my answer I said I assumed you wanted the mathematical result, which requires factoring to see if polynomials have common factors, and this then leaves the simplification problem. There is also the issue of how hard you want to look for common factors. But if you only want it to cancel preexisting factors, then it would be done in a different way that might miss potential mathematical common factors.

@epostma I certainly didn't think of that. In my purist view, that is not a unit - certainly my graduate students would never be allowed to write that in a paper, though I see a lot of students find it useful in their assignments.

@acer The math is OK, but the way it appears is certainly ugly. Even in my example, the cubic is messed up in a way that is no longer the conventional way to write it. I originally thought of of using factors on each of the numerators and denominators to try to keep things together, but that does not solve the expansion issue either. (I guess I long ago gave up using Maple to put things in a specifc nice form, because for some example it will fail.)

Edit: It seems that the way that factor normalizes as part of its operation is the central problem, which means it is hard to get around. factor((x+y)^3+x); fails to factor but leaves the answer in expanded form. So it could be improved by returning the original if it failed to factor, but that is probably a lot of additional cosmetic code.

@C_R Thanks - you provided a significant shortening with the nested subsindets, so a definite joint effort.

@C_R  The code could be incorporated into simplify using

`simplify/Unit` := ex -> subsindets(ex, specfunc(Units:-Unit), x -> subsindets(op(x), name, Unit))

and used as simplify(ex,Unit) (the normal seems to be automatic in this case.)

Expansion_of_compound_units_shortest.mw