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

@tomleslie Both sqrt(a) and a^b are uniquely defined and positive for any a > 0, b::real. The "cannot determine" issue isn't due to ambiguity; it's due to Maple avoiding floating-point evaluation of symbolic constant expressions unless the user forces it to with is or by some other means.

@bmartin Given that the OP wishes to access the posted procedure from a library, the with command (used at the top level) wouldn't be appropriate, and it doesn't work correctly if used within the procedure (as explained at ?with).

@nryq I put my code for the plot of #2 under my earlier Answer.

@nryq To plot #2, the sphere-capped cone, you need to know that the sphere's center is [0, 0, 1/2] and its radius is 1/2. This Answer doesn't cover the steps (completeing the square) to do that. It's just simple algebra that I find easier to do in my head than with Maple.

Cone:= z=sqrt(x^2+y^2):
Sphere:= z=x^2+y^2+z^2:
simplify(
    changecoords([Cone, Sphere], [x,y,z], cylindrical, [rho, phi, z]),
    symbolic
);
                    [                2    2]
                    [z = rho, z = rho  + z ]
plots:-display(
    plot3d~(
        [[z, phi, z], [sqrt(z-z^2), phi, z]],
        z=~ [0..1/2, 1/2..1], phi= -Pi..Pi, coords= cylindrical,
        scaling= constrained
    )
);

MaplePrimes isn't letting me copy-and-paste plots anymore! You'll need to run the above to see the plot. It looks like a squat ice cream cone.

@dharr The significant difference, which is obscured by your use of 2D Input, is that you've used t__1 instead of the OP's t[1].

@acer Thanks for finding the error. It's now fixed in the Answer above.