Maybe Robert Israel's sniffmem routine can help.

You could try these,

evalf( Int( a*r^2*'Phi'(r,10)^2 + 2*b*r*'Phi'(r,10)^3, r=0..infinity ) );

evalf( Int( r -> a*r^2*Phi(r,10)^2 + 2*b*r*Phi(r,10)^3, 0..infinity ) );

See the evalf/Int help page for additional options for numerical quadrature (numeric integration).

Could you add another equation to your original pde system, introducing a new term, the derivative (w.r.t r) of which is Phi(r,10)? ...or would it have to be Phi(r,t)? With the right BC, would solving that enhanced pde system provide a definite integral for Phi(r,10) w.r.t r?

> with(Statistics):

> Phiinverse := (x) -> Quantile(Normal(0, 1), x, numeric) :

> Phiinverse(1-1/2^24);

                          5.294718326

Note the accuracy, however,

> Digits:=16:

> evalf[10](Phiinverse(1-1/2^24));

                          5.294704085

> simplify(combine(Z,power));

                               0

Replace the square brackets inside the construction of newtonleast with round brackets.

Or, if you use insist on using square brackets there, then the angle-bracket constructor interprets it as a 1-by-2 Matrix and you would have to index its entries as newtonleast[1,1] and newtonleast[1,2].

In Maple, round brackets are operational delimiters and square brackets are list constructors.

You should change your code throughout. The procedure `IntermediateCalc` is extremely inefficient, especially in its heavy use of `unapply`. Use expressions when possible instead of all those procedures (almost all of which are not necessary).

If you get rid of the globals, and the lexical scoping, and type the parameters of the final procedure (which is the only one you really need to actually be a procedure) then you should be able to run it under evalhf or to compile it.

Instead of trying to run the Compiler to magically make it faster, you should first fix all the coding aspects that make it extra slow to begin with.

You should declare the arguments to your procedure with typechecking. You might want to change the type checks for `a` and `b` below from float to integer (and then pass them as such, accordingly, of course).

restart:
ggg_expression:=(2*(-a+6*y^2*a^2+2*b*y^2*a^2+9*b*y^3*a+3*b*y^2*a
+b*y^6*a^2+4*b*y^5*a^2+5*b*y^4*a^2+2*b*y^5*a+8*b*y^4*a+3*b*y^3*a^2
+4*b*y^2-y^3*a-3*y^2*a+y^6*a^3+5*y^5*a^3+8*y^4*a^3+4*y^3*a^3+y^5*a^2
+5*y^4*a^2+9*y^3*a^2+b*y^4+4*b*y^3+y*a^2+2*a*b*y+y^3*a*(y+1)^(-2*b)
+a*(y+1)^(-2*b)+3*y^2*a*(y+1)^(-2*b)+3*y*a*(y+1)^(-2*b)
-2*y^2*a^2*(y+1)^(-2*b)-y*a^2*(y+1)^(-2*b)-3*y*a-y^3*a^2*(y+1)^(-2*b))
/(y^2*(a+b)*(y+1)*(y+2)^2*(y*a+1)^2)+(-2-y^2+2*y^2*a^2-2*y+4*b*y^3*a
+8*b*y^2*a+2*b*y^4*a^2+4*b*y^3*a^2+2*b*y^2+4*b*y+2*y^2*a+y^4*a^2
+4*y^3*a^2+(y+1)^(-2*b)-2*y^2*a*(y+1)^(-2*b)-2*y*a*(y+1)^(-2*b)
+y^2*a^2*(y+1)^(-2*b)+(y+1)^(-2*b)*y^2+2*(y+1)^(-2*b)*y+(y+1)^(-2*b-2)*y^2
+2*(y+1)^(-2*b-2)*y+(y+1)^(-2*b-2)*y^4*a^2+2*(y+1)^(-2*b-2)*y^3*a
+2*(y+1)^(-2*b-2)*y^3*a^2+4*(y+1)^(-2*b-2)*y^2*a+(y+1)^(-2*b-2)*y^2*a^2
+2*(y+1)^(-2*b-2)*y*a+(y+1)^(-2*b-2))/(y^2*(a+b)*(y+2)^2*(y*a+1)^2))
/((-a+6*y^2*a^2+2*b*y^2*a^2+9*b*y^3*a+3*b*y^2*a+b*y^6*a^2+4*b*y^5*a^2
+5*b*y^4*a^2+2*b*y^5*a+8*b*y^4*a+3*b*y^3*a^2+4*b*y^2-y^3*a-3*y^2*a
+y^6*a^3+5*y^5*a^3+8*y^4*a^3+4*y^3*a^3+y^5*a^2+5*y^4*a^2+9*y^3*a^2
+b*y^4+4*b*y^3+y*a^2+2*a*b*y+y^3*a*(y+1)^(-2*b)+a*(y+1)^(-2*b)
+3*y^2*a*(y+1)^(-2*b)+3*y*a*(y+1)^(-2*b)-2*y^2*a^2*(y+1)^(-2*b)
-y*a^2*(y+1)^(-2*b)-3*y*a-y^3*a^2*(y+1)^(-2*b))^2
/(y^4*(a+b)^2*(y+1)^2*(y+2)^4*(y*a+1)^4)-(1/4)*(-2-y^2+2*y^2*a^2-2*y
+4*b*y^3*a+8*b*y^2*a+2*b*y^4*a^2+4*b*y^3*a^2+2*b*y^2+4*b*y+2*y^2*a
+y^4*a^2+4*y^3*a^2+(y+1)^(-2*b)-2*y^2*a*(y+1)^(-2*b)-2*y*a*(y+1)^(-2*b)
+y^2*a^2*(y+1)^(-2*b)+(y+1)^(-2*b)*y^2+2*(y+1)^(-2*b)*y+(y+1)^(-2*b-2)*y^2
+2*(y+1)^(-2*b-2)*y+(y+1)^(-2*b-2)*y^4*a^2+2*(y+1)^(-2*b-2)*y^3*a
+2*(y+1)^(-2*b-2)*y^3*a^2+4*(y+1)^(-2*b-2)*y^2*a+(y+1)^(-2*b-2)*y^2*a^2
+2*(y+1)^(-2*b-2)*y*a+(y+1)^(-2*b-2))^2/(y^4*(a+b)^2*(y+2)^4*(y*a+1)^4)):

ggg:=unapply(ggg_expression,[a::float,b::float,y::float]):

cggg:=Compiler:-Compile(ggg):

ggg(1.0,1.0,3.2);

                          9.924398803

cggg(1.0,1.0,3.2);

                      9.92439880079097314

N:=10000:

V1:=CodeTools:-Usage( Vector(N,(i)->ggg(1.0,1.0,2.0+1.0/i),datatype=float[8]) ):
memory used=0.88GiB, alloc change=39.37MiB, cpu time=20.44s, real time=20.62s

V2:=CodeTools:-Usage( Vector(N,(i)->evalhf(ggg(1.0,1.0,2.0+1.0/i)),datatype=float[8]) ):
memory used=1.90MiB, alloc change=0 bytes, cpu time=891.00ms, real time=899.00ms

V3:=CodeTools:-Usage( Vector(N,(i)->cggg(1.0,1.0,2.0+1.0/i),datatype=float[8]) ):
memory used=4.35MiB, alloc change=0 bytes, cpu time=187.00ms, real time=205.00ms

LinearAlgebra:-Norm(V1-V2), LinearAlgebra:-Norm(V2-V3);

                            -8                        -10
      3.23561790693815965 10  , 9.92976367797382409 10   

That piecewise is getting in the way. But if you take y strictly greater than zero, then you can get rid of the wrapping piecewise and make at least some kind of progress with series approximation.

expr:=MainF(y,a,b):

MultiSeries:-series(simplify(expr),y=0,7) assuming y>0:
R1:=simplify(convert(%,polynom)) assuming y>0:

simplify(expr) assuming y>0:
R2:=simplify(convert(series(%,y=0,7),polynom)) assuming y>0:

plot(eval([expr,R1,R2],[a=1,b=1]),y=0..0.4,legend=["MainF","R1","R2"]);

What form did you want, for this second example?

restart:

JA1S := (2*sqrt(2*y+3)*y+3*sqrt(2*y+3)-3*sqrt(3))/((2*y+3)^(3/2)*y):

sol := expand(frontend(simplify,[JA1S,size]));

                                 (1/2)    
                      1       3 3         
                      - - ----------------
                      y            (3/2)  
                          (2 y + 3)      y

simplify(JA1S-sol);

                               0

Applying evalf to the results from allvalues will actually cause Maple to unnecessarily duplicate effort (in the form of repeated calls to fsolve that re-do a bit of the earlier computational effort).

Basically, allvalues is making Maple internally use fsolve (with the avoid option) to get the floating-point "identifier" for the three distinct RootOfs. And then hitting those with evalf causes fsolve to be tasked once more (but this time with initial guesses of the identifiers). The second set of calls to fsolve (with RootOf identifiers as initial guesses) can be avoid by simply calling the first set of fsolve attempts manually.

You can see all those fsolve calls like so:

restart;
Digits:=15;
-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772;
E:=convert(%, rational);

trace(fsolve);
RootOf(E,x); [allvalues(%)]; evalf(%);

One way to do it manually, and which is the essence of what allvalue did above:

restart:
expr:=-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772:

S1:=fsolve(expr,x);

                          2.001200287

S2:=fsolve(expr,x,avoid={x=S1});

                          1.413550865

fsolve(expr,x,avoid={x=S1,x=S2}); # no solution, with unevaluated call as return value
 
      /                                  3  3                   
fsolve|-2.711505682 x - 7.65 ln(3 - x) - - x  + 8.422482772, x, 
      \                                  8                      

                                            \
  avoid = {x = 1.413550865, x = 2.001200287}|
                                            /

The above manual way (with the same sorts of calls to fsolve) can also be accomplished using Student:-Calculus1:-Roots.

restart:
expr:=-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772:

#trace(fsolve);  # uncomment to see the calls

Student:-Calculus1:-Roots(expr);

                   [1.413550865, 2.001200287]

Another quite different approach is to use RootFinding:-NextZero, although using it requires choosing a left-most point of the search.

restart:
expr:=-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772:
f:=unapply(expr,x):

S1:=RootFinding:-NextZero(f,0);

                          1.413550865

S2:=RootFinding:-NextZero(f,S1);

                          2.001200287

RootFinding:-NextZero(f,S2); # no solution

                              FAIL

Use fsolve intsead of solve here.

> restart:

> eq:=-2.711505682*x-7.65*ln(3-x)-3/8*x^3+8.422482772=0:

> fsolve(eq,x);

                          2.001200287

Raise Digits beforehand, for more digits in the result. After Digits:=90 all returned digits appear still to be correct for this example.

It should be cos(x)^2 not cos^2x

Your mistake is to not notice that the two results are equivalent. Hit your "correct" result with the expand command, or simplify their difference, etc.

I'll take the origin as the fourth (and common) point which is required for the geom3d specification of that object.

restart:

with(plots):
with(geom3d):

point(A,2,-3,1), point(B,-4,3,3), point(C,2,4,6), point(Origin,0,0,0):

dsegment(d1,[Origin,A]), dsegment(d2,[Origin,B]), dsegment(d3,[Origin,C]):

volume(parallelepiped(pp,[d1,d2,d3]));

                                      100

display( draw(pp,color=cyan),
         display(draw(Origin,symbolsize=25),color=red),
         display(draw([A,B,C],symbolsize=25),color=black),
              axes=box, scaling=constrained,
              view=[-4..4,-3..7,0..10],
              labels=[x,y,z] );

It is a bug, and it is present in 15.01.

Here is a way to see valid solutions for this example (in 15.01, but not tried in 13), by splitting across what seem to be relevant conditions.

restart:

myDeterminant := (1/2)*l[4]*l[1]*(l[4]*cos(x[2]-2*x[3])
                 +l[4]*cos(x[2])-l[1]*sin(2*x[2]-x[3])-l[1]*sin(x[3])):

simplify(myDeterminant,[l[1]=0]);
                               0

simplify(myDeterminant,[l[4]=0]);
                               0

S1:=solve({myDeterminant = 0,cos(x[2])<>0,sin(x[2])<>0,l[1]<>0,l[4]<>0}, AllSolutions):

S1:=subsindets(subsindets([S1],l[1]=l[1],t->NULL),l[4]=l[4],t->NULL);
     [ /                           1              \   
     [{ x[2] = x[2], x[3] = x[2] + - Pi + 2 Pi _Z1 }, 
     [ \                           2              /   

        /                           1              \   
       { x[2] = x[2], x[3] = x[2] - - Pi + 2 Pi _Z2 }, 
        \                           2              /   

        /       l[1] sin(x[2])                          \ ]
       { l[4] = --------------, x[2] = x[2], x[3] = x[3] }]
        \         cos(x[3])                             / ]

simplify(expand(simplify(eval~(myDeterminant,S1))));
                           [0, 0, 0]
S2:=solve({myDeterminant = 0,cos(x[2])<>0,sin(x[2])=0,l[1]<>0,l[4]<>0}, AllSolutions):

S2:=subsindets(subsindets([S2],l[1]=l[1],t->NULL),l[4]=l[4],t->NULL);
     [ /                          1              \   
     [{ x[2] = 2 Pi _Z5, x[3] = - - Pi + 2 Pi _Z7 }, 
     [ \                          2              /   

        /                        1              \   
       { x[2] = 2 Pi _Z3, x[3] = - Pi + 2 Pi _Z8 }, 
        \                        2              /   

        /                               1              \   
       { x[2] = Pi + 2 Pi _Z6, x[3] = - - Pi + 2 Pi _Z9 }, 
        \                               2              /   

        /                             1               \ ]
       { x[2] = Pi + 2 Pi _Z4, x[3] = - Pi + 2 Pi _Z10 }]
        \                             2               / ]

simplify(expand(simplify(eval~(myDeterminant,S2))));
                          [0, 0, 0, 0]

S3:=solve({myDeterminant = 0,cos(x[2])=0,sin(x[2])<>0,l[1]<>0,l[4]<>0}, AllSolutions):

S3:=subsindets(subsindets([S3],l[1]=l[1],t->NULL),l[4]=l[4],t->NULL);
[ /         1                                      \   
[{ x[2] = - - Pi + 2 Pi _Z14, x[3] = Pi + 2 Pi _Z17 }, 
[ \         2                                      /   

   /       1                                      \   
  { x[2] = - Pi + 2 Pi _Z13, x[3] = Pi + 2 Pi _Z18 }, 
   \       2                                      /   

   /         1                                 \   
  { x[2] = - - Pi + 2 Pi _Z12, x[3] = 2 Pi _Z19 }, 
   \         2                                 /   

   /       1                                 \   
  { x[2] = - Pi + 2 Pi _Z11, x[3] = 2 Pi _Z20 }, 
   \       2                                 /   

   /                              1                            \   
  { l[1] = l[4] cos(x[3]), x[2] = - Pi + 2 Pi _Z15, x[3] = x[3] }, 
   \                              2                            /   

   /                                 1                            
  { l[1] = -l[4] cos(x[3]), x[2] = - - Pi + 2 Pi _Z16, x[3] = x[3]
   \                                 2                            

  \ ]
   }]
  / ]

simplify(expand(simplify(eval~(myDeterminant,S3))));
                       [0, 0, 0, 0, 0, 0]

Less complete is to do it without AllSolutions. Eg.

restart:
myDeterminant := (1/2)*l[4]*l[1]*(l[4]*cos(x[2]-2*x[3])
                 +l[4]*cos(x[2])-l[1]*sin(2*x[2]-x[3])-l[1]*sin(x[3])):
solutions := solve({myDeterminant = 0,sin(x[2])<>0,cos(x[2])<>0});
nops([solutions]);
simplify(expand(simplify(eval~(myDeterminant,[solutions]))));
solutions2:= solve({myDeterminant = 0,cos(x[2])=0,sin(x[2])<>0});
nops([solutions2]);
{solutions2} minus {solutions};
nops(%);
simplify(expand(simplify(eval~(myDeterminant,[solutions2]))));
nops(%);
solutions3:= solve({myDeterminant = 0,cos(x[2])<>0,sin(x[2])=0});
nops([solutions3]);
#{solutions3} minus {solutions};
#nops(%);
simplify(expand(simplify(eval~(myDeterminant,[solutions3]))));
nops(%);
#{solutions} union {solutions2} union {solutions3};
#nops(%);