@dharr Thanks, I didn't know that convert(..., FormalPowerSeries) allowed an expansion point as an optional 3rd argument. 

@The function Yeah, "Pochhammer" definitely sounds like a heavy-metal band's name. But it's a fairly simple function similar to factorial.

@salauayobami Do you mean that you'd like to see a shooting method solution using RKF45 as the underlying IVP method? Or do you mean that you'd like some way of verifying that FDM is being used? 

@janhardo There's an entire textbook on Maple programming built right into Maple's help system. Just search "ProgrammingGuide" in the help browser. There are 20 such help "pages" (each chapter length). It includes much information on procedures and modules.

This is what I meant by the "bad programming practice" of returning values through the parameters, which is used by the procedures maxmin and optimize that you showed. To give a simpler example, let's say that I want a procedure that squares its input. Here is the bad way:

BadSqr:= proc(x, sqr::evaln)
    sqr:= x^2;
    return
end proc
:
BadSqr(3, s); #Executing this line produces no output.
s;
                               9

And here are two normal and acceptable ways of doing it:

BetterSqr:= proc(x)
    return x^2
end proc
:
BetterSqr(3);
                               9

Sqr:= x-> x^2:
Sqr(3);
                               9

The way that I call the bad way was considered bad long before 2002 (when that book was written).

@janhardo I don't know what you mean by this:

• I see that the input parameters from the second procedure as local variables are declared in the main procedure.
  Indeed the sub procedure is treated like a local variable.

Tom's procedure does not use maxmin at all, neither as an external procedure nor as a subprocedure of optimize. His local variables maxv and maxv could've been named anything. There is no connection to the procedure maxmin that you posted.

@salauayobami Maple's dsolve's numeric BVP solvers are all FDM. RKF45 is a method for IVPs. Occasionally an IVP solver (could be any) is used iteratively as part of a BVP solver in something called the "shooting method". Those are coded ad hoc by users.

Good job, @janhardo ! Vote up. Glad to see that you've learned some modern Maple coding.

@brian bovril There is nothing wrong with Tom Leslie's Answer. However, I'd like to point out that my Answer, while brief, both corrects the one and only error in your worksheet and gives a bit of background information regarding the source of the error, which is understandable given the likelihood of mixing up G and g. So don't you think that it deserves a vote up also?

If my Answer had been posted after Tom's, then indeed it would seem somewhat superfluous. But it wasn't; I posted it 6 hours before Tom's, less than an hour after you posted the Question.

@dharr Your correction of the boundary condition D(f)(0)^2 to (D@@2)(f)(0) is sufficient to correct the "matrix is singular" error, even after the multiplication after delta is used.

@dharr In ODE1, delta appears as a function symbol, but it was intended to be a coefficient. It needs to be followed by a multiplication.

This will lead to a "matrix is singular" error from dsolve. My first try for dealing with that would be to pertube the parameters a little.

I'm not sure, but I think that "too many levels of recursion" is a form of stack overflow, one of the stated types of untrappable error.

@janhardo There are two differences between those two versions of Gradient:

1. Return type: The 1st returns a procedure that returns a list; the 2nd returns a list of procedures. The 2nd is slightly more convenient when calling fieldplot.

2. Arity: The 1st is limited to functions of two variabes; the 2nd works for any number of variables and defaults to 2. (The number of variables of a function is called its arity.)

@Anthrazit Are you saying that the version used in current Maple (Maple 2021) was considered obsolete 6 years ago? 

I updated (below) the code given in my Answer above. I added the following:

1. A paragraph comment at the top of each procedure.

2. Filtering out of nonreals as suggested by Acer.

3. Critical points are color-coded, and the same color used in a legend in the table (via MathML).

4. The LCHuv color space is used to achieve maximal perceptible difference among the selected colors.

5. The surface plot is (by default) given in surfacecontour form.

6. A 2D fieldplot of the gradient is embedded as the "floor" of the 3D plot.

7. The output of the final procedure can be the plot, the table, or both (the default). If both are wanted, a Record is returned.

restart
:
#Returns a list of procedures that return the 1st partial derivatives #of f.

Gradient:= proc(f, {arity::posint:= 2})
local x, i, V:= seq(x[i], i= 1..arity); 
    map(unapply, diff~(f(V), [V]), [V])
end proc
:
#Returns a procedure that returns a symmetric matrix of the 
#2nd partial derivatives of f. The matrix is declared symmetric so that
#its eigenvalues will be computed as strictly real.

Hessian:= proc(f, {arity::posint:= 2})
local x, i, V:= seq(x[i], i= 1..arity);
    unapply(
        Matrix(arity, (i,j)-> D[i,j](f)(V), 'shape'= 'symmetric'), [V]
    )
end proc
:
#Convert a color from plot form to *ML form.

PlotColorToML:= ColorTools:-RGB24ToHex@ColorTools:-PlotColorToRGB24
:
#Converts a string to a symbol in MathML form that'll display in an
#upright font in black, possibly with a background color.

Mn:= proc(s, C::specfunc(COLOR):= ())
local 
    c:= if C=() then 
        else `, mathbackground= "` || (PlotColorToML(C)) || `"`
        fi
;   
    `#mn("` || s || `", mathcolor= "#000000"` || c || `)`
end proc
:
#2nd Derivative Test: Use H (the Hessian matrix) to classify a critical
#point as minimum, maximum, or saddle point. This works for any 
#number of dimensions, whereas the well-known formula that uses the
#determinant of the Hessian only works for 2-variable functions. 

Classify:= proc(H::procedure, P::list(realcons))
local E:= LinearAlgebra:-Eigenvalues(evalf(H(P[])));
    if member(0., E) then "inconclusive" 
    elif max(E) < 0 then "max"
    elif min(E) > 0 then "min"
    else "saddle"
    fi
end proc
:
#Find the zeroes of the gradient. No attempt is made to handle 
#discontinuities, nondifferentiability, or parameterized solutions 
#returned by `solve` (_Zn, _Bn, _NNn, etc.).

CritPts:= proc(f)
local x, y, V:= (x,y);
    select(
        p-> Im~({p[]}) = {0.}, #Filter out nonreals.
        (`simplify/zero`@fnormal@map2)(
            eval,
            [V, f(V)], 
            [evalf@solve](
                {Gradient(f)(V)[]}, {V}, 'explicit', 'allsolutions'
            )
        ) 
    )
end proc
:    
#Create some parameter types:

AddType:= (t::name, p::{type, procedure})->
    if not TypeTools:-Exists(t) then TypeTools:-AddType(t,p) fi
:
#half-open unit interval:
AddType('`[0,1)`', And(nonnegative, negative &under (`-`, 1)))
:
#list of options:
AddType('Opts', list(name= anything))
:
#Given a list L of reals, return a range that contains L plus margins
#on each side such that (total margin width)/(range width) = `margin`.

ViewRange:= proc(L::{set,list}(realcons), margin::`[0,1)`:= 0.3) 
local M,m,d;
    (M,m):= (max,min)(L); 
    d:= margin*(M-m)/2/(1-margin);
    (m-d)..(M+d)
end proc
:
#Choose a bright color scaled for maximal perceptible separation and
#return in plot-color form. The magic numbers 68 and 122 come from the
#average values of the L and C coordinates for 
#    ColorTools:-Convert([x, 1., 1.], "HSV", "LCHuv"), x= 0..1.
#For more info, see the Wikipedia article "HCL color space"
#https://en.wikipedia.org/wiki/HCL_color_space
#
#Some valid LCHuv combinations aren't representable in RGB. For those
#cases, the adjustment max~(0, min~(1, [R,G,B])) is used. 

ChooseColor:= proc(C::`[0,1)`:= rand(0. .. 1.)())
    'COLOR'(
        'RGB',
        max~(0, min~(1, 
            ColorTools:-Convert([68., 122., 360.*C], "LCHuv", "RGB")
        ))[]
    )
end proc
:
#Create a plot of the surface with color-coded and indexed critical
#points, a field plot of the gradient, and a data table listing and
#classifying the critical points with a legend keyed to the color
#coding and indexing in the plot.

PlotAndTable:= proc(
    f, 
    V::list(name):= [x,y,z],
    {
        margins::And(list(`[0,1)`), 3 &under nops):= [0.4, 0.4, 0.2],  
        pointopts::Opts:= [
            'symbol'= 'solidsphere', 'symbolsize'= 16
        ],
        funcopts::Opts:= [
            'style'= 'surfacecontour', 'transparency'= 0.5
        ],
        textopts::Opts:= [
            'align'= {'ABOVE', 'RIGHT'}, 
            'font'= ['courier', 'bold', 24], 'color'= 'black'
        ],
        fieldopts::Opts:= ['arrows'= 'SLIM', 'fieldstrength'= 'log'],
        output::identical(plot, table, record):= record
    }
)
uses P= plots;
local 
    C:= CritPts(f), n:= nops(C), G:= Gradient(f), 
    VR:= V=~ViewRange~(map2(index~, C, [1,2,3]), margins), 
    k, L:= <seq(1..nops(C))>, 
    Colors:= ChooseColor~([seq](k/n, k= 0..n-1)),
    R:= Record("plot"= (), "table"= ())
;
    R:-plot:= P:-display(
        plot3d(f(V[]), VR[1..2][], funcopts[]),
        P:-pointplot3d(C, color= Colors, pointopts[]),
        P:-textplot3d(
            {seq}([C[k][], cat(" ",k)], k= L), textopts[]
        ),
        #This is a trick to embed a 2D plot as the bottom plane in a 3D
        #plot:
        (PLOT3D@subsindets)(
            op(1, P:-fieldplot(G, rhs~(VR[1..2])[], fieldopts[])),
            [numeric $ 2], (L,z)-> [L[], z], op([2,1], VR[3])
        ),
        'view'= rhs~(VR), 'labels'= V,
        _rest
    );
    R:-table:= < 
        <'Legend' | V[] | 'Type'>, 
        < 
            Mn~(cat~(" ", L, " "), Colors) | 
            Matrix(C) | 
            map2(Mn@Classify, Hessian(f), <C>)
        >
    >;
    `if`(output=record, R, R[output]) 
end proc
:

Usage:

R:= PlotAndTable((x,y)-> x^4 - 3*x^2 - 2*y^3 + 3*y + x*y/2);

R:-plot;

MaplePrimes won't let me copy-and-paste plots! See the worksheet.

R:-table;

Of course, the table columns are properly aligned in the worksheet!

Download ExtremaAndSaddles.mw

@maxbanados I substantially updated my Answer above at the same time as you were asking your followup Question about integrals. So, that's worth another reading.

For your integrals, the operation that you specified only makes mathematical sense to me for definite integrals, so that's what I coded it for. However, it could easily be changed to handle indefinite integrals or both.

Here's pseudocode of my algorithm. Almost any aspect of this can be easily changed if you want.

Algorithm CollapseInts:
Input: An expression e possibly containing integrals to convert; function names f and g.

1. Select from e all definite integrals whose integrand is f(...) (rather than simply contains f(...)). The fs may have any number of arguments, which may be any expressions, not necessarily simple names.

2. FOR every integral J in the set selected in 1 DO

1. Let f1 be the integrand and d the differential variable.
2. IF f1 contains d as an argument in any single position AND
       the remaining arguments of f1 have no dependence on d
   THEN replace J in e with g applied to the arguments of f1 with d removed.

And here's the Maple code:

CollapseInts:= proc(e, f::name, g::name)
    subsindets(
        e,
        And(specfunc({Int, int}), anyfunc(specfunc(f), name=range)),
        proc(J)
        local 
            d:= op([2,1], J), #differential variable
            f1:= op(1,J), #f(...)
            p #position of d in f1
        ;
            if member(d, f1, p) and 
                not has([op](f1)[[1..p-1, p+1..-1]], d)
            then g(op(subs(d= NULL, f1)))
            else J
            fi
        end proc
    )
end proc
:
#Example usage:
e:= Int(f(x,y), y= a..b)+Int(1+f(x,y), y= a..b)+Int(f(1+y,y), y= a..b)+Int(f(y,1+x), y= a..b):
CollapseInts(e, f, g);
       /  /b               \   /  /b               \           
       | |                 |   | |                 |           
g(x) + | |   1 + f(x, y) dy| + | |   f(1 + y, y) dy| + g(1 + x)
       | |                 |   | |                 |           
       \/a                 /   \/a                 /           

Using 1D input (aka Maple Input) and a fairly recent version of Maple, the above can be shortened to

CollapseInts:= (e, f::name, g::name)->
    subsindets(
        e,
        And(specfunc({Int, int}), anyfunc(specfunc(f), name=range)),
        J-> 
        local d, f1, p;            
            if member((d:= op([2,1], J)), (f1:= op(1,J)), p) and
                not has((f1:= subsop(0= g, p= (), f1)), d)
            then f1
            else J
            fi
    ):