@jennierubyjane You can include fixedparameter in the eval call in the dsolve call by replacing Pr= L[k] with

{Pr= L[k], fixedparameter[]}

But you can't use Nb=0 because Nb is a denominator in one of the ODEs.

Here's an interesting fractalish plot:

A:= Iterate(x-> 3.7*x*(1-x), 0.23, 10^7): #using same Iterate defined above
T:= Statistics:-TallyInto(A, 0..1, bins= 1000):
PlotIt:= proc(T::list(range(realcons)= nonnegint))
local r:= op~(1,T), L:= lhs~(r), R:= rhs~(r), F:= op~(2,T);
    plot(
        `[]`~(`[]`~(L,0), `[]`~(L,F), `[]`~(R,F), `[]`~(R,0)),
        color= black, thickness= 0
    )
end proc
:
PlotIt(T);

@janhardo After differentiating the evaluated integral, you can get back to the original integrand like this:

J:= Int(exp(t*x), x);
J1:= subs(t= a+b*I, J);
evalc(J1);
value(%);
simplify(convert(diff(%, x), exp), {a+b*I= t});
                            exp(x t)

The order of applying evalc and value doesn't matter.

@janhardo When you use the notation x+I*y, you're likely assuming, perhaps subconsciously, that x and y are real; however, Maple doesn't automatically assume this. It assumes that all variables are complex unless told otherwise, and regardless of the presence of I. The command evalc will temporarily assume that variables are real so that expressions can be separated into their real and imaginary parts. So, compare Re(x+I*y) with evalc(Re(x+I*y)).

It's not needed to explicitly have an I for complex operations to be used. They're simply used by default. For example, the derivative of sin(x) is cos(x) regardless of whether x is real. Unless told otherwise Maple is not assuming that it's real.

@Pepini The code that you showed doesn't actually use TallyInto; it uses Histogram instead. That's a fine alternative if you don't need programmatic access to the actual counts for some future purpose.

Here's complete code for your problem:

Iterate:= proc(f, x0::complexcons, N::posint)
local A:= Array(1..N, [x0], datatype= hfloat);
    evalhf(
        proc(f, A, N)
            for local k from 2 to N do
                A[k]:= f(A[k-1])
            od
        end proc
        (f, A, N)
    );
    A
end proc
:
A:= CodeTools:-Usage(Iterate(x-> 3.7*x*(1-x), 0.2, 10^7)):
memory used=76.30MiB, alloc change=76.30MiB, 
cpu time=2.83s, real time=2.85s, gc time=0ns

T:= Statistics:-TallyInto(A, 0..1, bins= 40);
PlotIt:= proc(T::list(range(realcons)= nonnegint))
local x;
    plot(
        [seq]([
            [op([1,1],x),0], [op([1,1],x),op(2,x)], 
            [op([1,2],x),op(2,x)], [op([1,2],x),0]
        ], x= T
        ), color= black, thickness= 0
    )
end proc
:
PlotIt(T);

It takes about 3 seconds to generate the 10 million points. The time for the tallying and plotting is too small to even measure.

@nm Please let me know if you find any cases where the identity method doesn't work, either false positive or false negative. I'd be surprised if there was a false positive (an incorrect p value returned), but not surprised at a false negative (no p found when the function is indeed isobaric). For any false negative, I wonder if any other general technique can find the p.

The documentation for solve(identity(...)) says that the syntax allows for the declaration of one free variable. I chose a to be the declared free variable, but x and y are also free. Thus, my confidence in this technique is not as high as it could be.  

Yes, the sentence that you mentioned should be "Define new dependent variable u(t)."

I'm using a QHD display (3200 x 1800) with the Zoom factor in the Windows settings set to the recommended value of 250%. This is independent of Maple's own Zoom factor.

@Christopher2222 Yes, when it happens for me, it affects all open documents and worksheets.

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