@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?

The functions of x are just dead weight in your computation because there are no differential operations with respect to x. You might as well just make your equation

fracdiff(u(t), t, 1/2) = sqrt(t)/GAMMA(3/2)

That'll make it easier to solve regardless of the method you use (Maple, other CAS, by hand).

Of course, although this equation is in a sense trivial, it's still beyond the documented capabilities of intsolve.

@Kitonum That's a good and interesting example. Note also that add('1,2,3') works, showing that the issue is specifically evaluation, not directly the presence of a variable.

@acer Thanks for the tip. I was refining the whole code anyway, and I had already added option 'real' to the solve command, although I'm not completely confident that that'd eliminate the nonreal problem. To be safer, I'll also filter the floats for nonreals. What do you think?

@janhardo Before I expend too much effort trying to explain the procedures, I need to know whether you understand the math behind this, So, please describe in your own words the method of classifying critical points of multivariate functions known as the "Second-Derivative Test".

Re finding good plot ranges: We have a list C of critical points, each with x, y, and z values. Let's say we want the x range. Extract the first coordinate from each point in C. Let m be the minimum of these, and let M be the maximum. We don't quite want to make the range just m to M; we want some margin on each side of that. So, we make the range (m - r*(M-m)) .. (M + r*(M-m)) where 2*r is the desired fraction of the overall width that should be margin.