[This is an Answer to your followup Question "DirectSearch skips a parameter".]

There are several issues:

1. Neither I nor anyone else here can duplicate what you're doing because your original Question doesn't show the value of dsys, your system of ODEs. So, what I've said so far has necessarily been based on guessing.

2. DirectSearch:-DataFit has 9 different ways of computing residuals, including your standard least squares. Perhaps one would suit your needs?

3. Your resid needs to politely ignore non-numeric input the same way your ffit2 does. Also, you should use add instead of sum (sum is intended for symbolic summation).

resid:= (eta_L, L0)->
    if [args]::list(numeric) then add((ffit2~(eta_L, L0, X0) - Y0)^~2)
    else 'procname'(args)
    fi
:

4. This was probably just a transcription error as you were moving text to MaplePrimes, but note that you used resid2 in the code below your definition of resid.

A radical conversion is not always possible (usually not possible for irreducibles with degree > 4), but viewing the RootOf form of the factors is also quite unwieldy. Here's a way that I think helps one to see the relationships between the irrational roots:

restart:
p:= x^4 - x^2 + x - 1:

(* This procedure is only to compactify the display of RootOf
expressions. It has no effect computationally; in particular, further
processing of output displayed by this should ignore the radical symbol
(`√`) and Z and treat them as the usual RootOf and _Z.
*) 
`print/RootOf`:= p-> local Z; `√`(subs(_Z=Z, p), _rest):

(* This command--equivalent to the 1st step of Kitonum's Answer--shows
the factorization with nested RootOfs, from which it's still unwieldy 
to understand the relationships between the irrational roots.
*)
evala(AFactor(p)); #equivalent to PolynomialTools:-Split
                                     /           /  
          /           / 3    2    \\ |           | 2
  (x - 1) \x - √\Z  + Z  + 1// \x - √\Z 

                                                          2
       /           / 3    2    \\            / 3    2    \ 
     + \1 + √\Z  + Z  + 1// Z + √\Z  + Z  + 1/ 

                           \\ /                               
              / 3    2    \|| |               / 3    2    \   
     + √\Z  + Z  + 1/// \x + 1 + √\Z  + Z  + 1/ + 

           /                                 
           | 2   /           / 3    2    \\  
    √\Z  + \1 + √\Z  + Z  + 1// Z

                           2                       \\
              / 3    2    \           / 3    2    \||
     + √\Z  + Z  + 1/  + √\Z  + Z  + 1///

#Of course, the above looks better in Maple's 2D output, which I can't 
#display here. `√` (HTML code) will be display as a single-character-width root
#symbol.

(* This technique recursively nests and 'veils' the RootOfs, which I
think makes it easier to see the relationships. Here, nu could be any
unassigned name.
*)
LE:= LargeExpressions:
evalindets(evala(AFactor(p)), RootOf, LE:-Veil[nu]), 
    <seq(nu[k]=LE:-Unveil[nu](nu[k]), k= 1..LE:-LastUsed[nu])>
(x - 1)*(x - nu[1])*(x - nu[2])*(x + 1 + nu[1] + nu[2]), 
nu[1] = `&radic;`(Z^3 + Z^2 + 1), 
nu[2] = `&radic;`(Z^2 + (1 + nu[1])*Z + nu[1]^2 + nu[1])

Hint: The slope of any line equals the trigonometric tangent of the angle that it makes with the positive x-axis. The slope of a tangent line to a curve equals the curve's derivative at the point of tangency.

op(1,A) will give you the _Z^2+1. I don't think that there's a way to get back the original variable (assuming there was one) other than direct substitution:

subs(_Z= x, op(1,A))

In your specific case, you could use op(A) instead of op(1,A), but it's not a good idea in general.

Like this:

P:= y-> [-y^2, y]:
dist:= (A::list, B::list)-> sqrt(add((A-~B)^~2)):
solve(diff(dist(P(y), [0,-3]), y), y, real);
                               -1
P(%);
                            [-1, -1]

Maple does have general implicit plotting commands for equations in 2 or 3 variables (their names are plots:-implicitplot and plots:-implicitplot3d). However, for curves that can be expressed as polynomials in 2 variables, there's another command that gives a much better plot: algcurves:-plot_real_curve.

#Express curve as an expression implicitly equated to 0 rather than as
#an equation:
C:= y*(y^2 - 1)*(y - 2) - x*(x - 1)*(x - 2):
algcurves:-plot_real_curve(C, x, y);

 

dydx:= -diff(C,x)/diff(C,y): #implicit derivative
#x-coordinates of horizontal tangent points, exactly:
eval~(x, {solve}({C, dydx}, {x,y}, explicit, real)); 
                  /    1  (1/2)      1  (1/2)\ 
                 { 1 - - 3     , 1 + - 3      }
                  \    3             3       / 

evalf(%); #approximations to above:
                  {0.4226497307, 1.577350269}

#Tangent line at P= [x0,y0]:
T:= P-> eval(y = eval(dydx, [x,y]=~ P)*(x - x0) + y0, [x0,y0]=~ P):
T([0,1]);
                           y = -x + 1
T([0,2]);
                              1      
                          y = - x + 2
                              3      

I'm not sure that this Answer will competely solve your problem; however, it's obvious that what you currently have (in your first Reply) can't possibly work because L0= 1e10 is not in the interval L0= 1.5e10..1.7e10.

Indeed, the work integral for any tank whose horizontal cross sections are circles can be done by a single half-line formula, which is given in procedure PumpWorkInt below. I did parts a, b, c, and I added a2 (upright cone) and c2 (cap hemisphere) to show that turning a and c over doesn't change their volume, but makes a huge difference in the amount of work.

restart:

PumpWorkInt:= proc(
    r::procedure, #radius= f(height) formula 
    hrange::range(algebraic) #e.g., h= 0..H    
)
local H:= rhs(hrange), h, W:= Int(rho*g*Pi*r(h)^2*(H-h), h= hrange);
    W = value(W) assuming H > 0
end proc
:
Volume:= proc(r::procedure, hrange::range(algebraic))
local H:= rhs(hrange), h, V:= Int(Pi*r(h)^2, h= hrange);
    V = value(V) assuming H > 0
end proc
:
Tanks:= ((Record@op)~@map2)(
    `~`[`=`], ["name",       "radius",       "hrange"],
    <["a1: inverted cone",   h-> R/H*h,          0..H],
     ["a2: upright cone",    h-> R/H*(H-h),      0..H],
     ["b: sphere",           h-> sqrt(H^2-h^2), -H..H],
     ["c1: bowl hemisphere", h-> sqrt(H^2-h^2), -H..0],
     ["c2: cap hemisphere",  h-> sqrt(H^2-h^2),  0..H]
    >
);

Vols:= <seq(
    <nprintf(`#mn("%s")`, T:-name) | Volume(T:-radius, T:-hrange)>,
    T= Tanks
)>;
               
CommonVals:= 
    g= 981/100*Unit(m/s^2), rho= 1000*Unit(kg/m^3), R= 1/2*Unit(m)
:
Work1:= <
    Vols[.., 1] | <seq(PumpWorkInt(T:-radius, T:-hrange), T= Tanks)>
>:    
Work2:= <
    Work1 |
    combine(
        eval~(
            rhs~(Work1[.., -1]), 
            [seq]([CommonVals, H= h*Unit(m)], h= (2,2,1,2,2))
        ),
        units
    )
>:
Work3:= <Work2 | evalf[3](convert~(Work2[.., -1], units, kJ))>;

 Download PumpWork.mw

I don't know whether you consider a typed command to be doing it "manually", or if that means via pen and paper to you. So, here's how to do it with a typed command: You say that you've "been given values of the quotient C__A0/C__AA". I will suppose that that given value is q. Then do

simplify(solve({-k*t = 1/C__A0 - 1/C__AA}, k), {C__A0/C__AA= q});
                        /      q - 1  \ 
                       { k = --------- }
                        \    C__AA t q/ 

I don't know how to do this (or much else) with menu-based commands. The presence of that side relation with q adds a level of complexity to menu-based iteractions.

You can create a custom type Monomial like this:

TypeTools:-RemoveType('Monomial'): #avoid silly warning
TypeTools:-AddType(
    'Monomial', #quotes protect against global reassignment of names
    e->
        local Var:= 'Y'[integer $ 2], Pow:= {Var, Var^rational};
        e::Pow or
        e::'specop'(
            {Pow, Not({constant, 'satisfies'(e-> hastype(e, Var))})},
            `*`
        ) 
)
:
#Usage:
type(Y[1,-3]^(-7/4)*sqrt(t+x), Monomial);
                              true

type(Y[4,3]^3*sin(Y[1,2]), Monomial);
                             false

To understand the above, you should read the help pages for the types integer, fraction, and constant, the type operator satisfies, and the commands AddType and hastype. The type `*` and the type operators { },  specop, and Not are documented (minimally) on the page ?type,structure.

The above will not work in 2D Input or in some older versions of Maple. If you're using either of those, let me know.

When entering an integral as a Maple command (rather than simply for display), do not use "d". The variable of integration is determined by the variable immediately after the comma. The pretty-printed output will still show a "d" followed by a variable.

Use this:

Exponent:= (x::algebraic)-> 
    `if`(x=0, undefined, `if`(x=1, 0, `if`(x::{`^`, specfunc(`%^`)}, op(2,x), 1)))
:

So, the result is undefined if the term is 0, it's 0 if the term is 1, it's 1 if the term is not 0 or 1 but has no explicit exponent, and it's the explicit exponent if one appears.

Note that non-numeric denominators always become multiplicands with their exponent negated in their internal representation.

There are two problems. The first is a syntax error in one of the boundary conditions. You've entered

varepsilon + delta*D@@2*(f)

You need to change that to

varepsilon + delta*(D@@2)(f)(LB)

The second problem is that you can't use infinity as the upper boundary UB; the current numeric BVP algorithms won't allow it. I changed UB to 5. I also changed Digits to 40 because a warning message advised raising it to 32. With these changes, I got through the first set of plots. Using different parameter values, it's very likely that you'll get to some boundary-layer effects, for which you'll get the error message "Newton iteration is not converging" or "initial Newton iteration is not converging". Fixing this problem (if it happens) is a much more subtle process. 

The Sieve of Eratosthenes is great for finding large sets of small primes, but it's useless (much too slow) for finding large primes. It seems like you're interested in finding large primes, such as those found by GIMPS. Below, I have implemented the Lucas-Lehmer test for Mersenne primes, which is the same algorithm used by GIMPS.

#Compute N mod (2^n-1) via bitwise operations. N must be < 4^n.
BitMod:= (N::posint, n::posint)->
    if ilog2(N+1) < n then N
    else `if`(ilog2(N) < n, 0, thisproc(add(Bits:-Split(N, n)), n)) 
    fi
:
#Check whether (2^p-1) is prime by Lucas-Lehmer algorithm.
IsMersennePrime:= proc(p::And(prime, Not(2)))
local s:= 4;
    to p-2 do s:= BitMod(s^2, p) - 2 od;
    evalb(s=0)
end proc:
IsMersennePrime(2):= true
:
#This Mersenne prime was discovered by computer in 1979:
CodeTools:-Usage(IsMersennePrime(44497));
memory used=2.65GiB, alloc change=0 bytes, 
cpu time=10.58s, real time=12.60s, gc time=6.48s

                              true
#Verify a negative case:
nextprime(44497);
                             44501

CodeTools:-Usage(IsMersennePrime(44501));
memory used=2.64GiB, alloc change=0 bytes,
cpu time=10.34s, real time=12.54s, gc time=6.39s

                             false
#Compare timings using default method not specific to Mersenne numbers:
restart:
CodeTools:-Usage(isprime(2^44497-1));
memory used=3.73GiB, alloc change=61.50MiB, 
cpu time=59.81s, real time=63.96s, gc time=1.92s

                              true

restart:
CodeTools:-Usage(isprime(2^44501-1));
memory used=2.78MiB, alloc change=0 bytes, 
cpu time=21.45s, real time=21.47s, gc time=0ns

                             false

To make your example work, all that you need to do is remove the export from in front of move in module dog. In other words, this works:

animal:= module()
option object;
export 
    data1,
    move::static:= proc(_self, $)
        print("In Animal class. moving ....")
    end proc
;
end module;
 
#create class/module which extends the above
dog:= module()
option object(animal);
    move::static:= proc(_self, $)
        print("In dog class. moving ....")
    end proc  
end module;

You can verify that the method overriding and polymorphism has been done properly.

To make that work in library code, the redefinition of move should be put in a ModuleLoad procedure, like this:

dog:= module()
option object(animal);
local 
    ModuleLoad:= proc()
        move::static:= proc(_self, $)
            print("In dog class. moving ....")
        end proc;
        return
    end proc
;
    ModuleLoad()  
end module;