@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
    ):

@janhardo Okay, I now understand your confusion regarding gradients.

Almost everything in calculus can be interpreted in a geometric/graphical way and an algebraic/symbolic way. The geometric way helps with "seeing the big picture"; the algebraic way is easier to compute with. Your understanding of gradients is geometric: The gradient is the vector that points in the direction of locally maximal steepest ascent, and it's perpendicular to the level curves. The algebraic interpretation is that the gradient is the vector (or list) of 1st partial derivatives, and it equals 0 at the critical points. Likewise, the Hessian is the matrix of all 2nd partial derivatives. The signs of its eigenvalues can be used to classify the critical points. In the case of 2-variable functions f(x,y), the determinant can be used instead of the eigenvalues.

@John2020 It's odd that this use of the documented option laplace doesn't work:

fracdiff(u(t), t, 1/2, method= laplace)

yet this does work:

inttrans:-laplace(fracdiff(u(t), t, 1/2), t, s))

And both statements remain true if you replace u(t) by u(x,t).

You want to change all F(f(x)) to H(h(x)) "where f(x) is really 'anything'". Like @mmcdara , I'm having trouble understanding that. We need to be able to extract the x from f(x) in order to place it in h(x). You can't extract x from simply anything at all.

@nm I only mention the following things for the sake of convenience. There isn't anything wrong with your Answer.

1. For several years now, specfunc(f) has been an acceptable abbreviation for specfunc(anything, f).

2. Both specfunc and anything are initially protected names. Thus, I feel less need to put them in quotes. However, many type names are not protected.

3. specfunc(f) matches functions named f with any number of arguments, including 0.

@nryq The part in red is the Maple code. I believe that this code will work in any Maple version since Maple 13, in any user interface or input mode. The part in blue is the output of the first section of code. There's no need to transcribe that.

Only the 2nd paragraph of the Maple code is needed to get the plot. The purpose of the 1st paragraph of code is simply to give you a sense of how I derived the parametric cylindrical-coordinate representations of the cone and sphere that I used, specifically

• Cone: [z, phi, z] (transformed from z = sqrt(x^2 + y^2))
• Sphere: [sqrt(z - z^2), phi, z] (transformed from z = x^2 + y^2 + z^2)

What is your mathematical familiarity with alternate coordinate systems such as polar, spherical, cylindrical?