The relevant equations are given on the web page "Effect of atmospheric drag on artificial satellite orbits". The 2nd equation (10.131) gives the impression that the atmospheric density is a simple function of altitude. That's only an approximation; the density is also affected by the presence of magnetic storms and how much sunshine the specific portion of atmosphere is receiving.

The plot below is almost identical to your model plot. The value of M is 1408, i.e., the number of values of r used, 3 <= r <= 4. These values are chosen by the adaptive routine of plot. The value of iter ranges from a low of 110 to a high of over 450,000, depending on r. It is chosen by my procedure for convergence of the Lyapunov limit to 5 decimal places. It takes about 16 seconds to produce the plot.

I used my procedure from 5 years ago that Acer posted a link to. It only required one minor modification due to the extremely strict syntax requirements of evalhf.

LyapunovExponent:= proc(
     F::procedure,
     x0::complexcons,
     {epsilon::positive:= 1e-7}, #absolute error tolerance
     {max_iters::posint:= 2^21}   #infinite loop prevention
)
option `Author: Carl J. Love, 2016-May-10`;
description
"Computes a numeric approximation to "
"lambda(x0) = limit(sum(ln(abs(f'(x[i]))), i= 0..n-1)/n, n= infinity), "
"where x[i] = f(x[i-1]), x[0]= x0. "
"See: https://en.wikipedia.org/wiki/Lyapunov_exponent"
;
local 
    #Trick evalhf into believing that F has no lexical variables:
    f:= parse(sprintf("%a", eval(F))), 
    f1:= D(f)
;
    evalhf(
        proc(f, f1, x0, epsilon, max_iters)
        local S:= 0, L:= epsilon, newL:= 0, x:= x0, n;
            for n to max_iters while abs(newL-L) >= epsilon do
                L:= newL; 
                newL:= (S:= S+ln(abs(f1(x))))/n; 
                x:= f(x)
            od;
            if n > max_iters then
                error
                    "didn't converge; abs. err. = %1; "
                    "increasing epsilon or max_iters might help.",
                     abs(newL-L)
            fi;
            userinfo(
                1, LyapunovExponent, "Converged; iterations =", n
            );
            newL
        end proc
        (f, f1, x0, epsilon, max_iters)
    )
end proc
:
f:= a-> x-> a*x*(1-x):
infolevel[LyapunovExponent]:= 1:
CodeTools:-Usage(
    plot(
        a-> LyapunovExponent(f(a), 0.2, epsilon= 1e-5), 3..4, 
        smartview= false
    )
);

MaplePrimes will not let me post the plot (I guess due to the fairly high number of points).

It can be done by

Ont_data:= rawdata[rawdata[prname] =~ "Ontario"];

See the help page ?DataFrame,indexing, particularly the "Boolean" sections.

The command is timelimit. You had timeout.

If you add option implicit to dsolve, this ode is solved instantly.

The number of boundary conditions required is the differential order + the number of symbolic parameters (those without numeric values). The differential order is the sum of the orders of the highest-order derivatives of each function; so, that's 7: 3 for f, 2 for theta, and 2 for phi. You gave a numeric value to the parameter Pr. The remaining symbolic parameters are Bi, Nb, Le, and Nt. You need to either give numeric values to 3 of these 4, or come up with 3 more boundary conditions.

Also, you need to change the square brackets to parentheses in the 2nd boundary condition in BC2. Square brackets can never be used in place of parentheses; they have their own meaning in Maple.

The command Statistics:-TallyInto does what you want.

The complex numbers are the default number system in Maple, and most (not all) commands (other than plots) that work with real-valued functions of real variables will also work with complex-valued functions of complex variables without needing any change in syntax or the addition of any options.

You wrote:

• It seems that there are a lot of standard calculus statements can be used by adding the word : complex 

I'm aware of limit and fsolve having a complex option. What else have you found?

Your expression is a polynomial. Here is a procedure to get all of its coefficients, including the constant term:

CoeffsOfIndexed:= proc(p, X::name)
local T, V:= indets(p, X[anything]), C:= [coeffs](p, V, 'T'), t;
    if p::linear(V) then
        table(sparse, [seq](`if`(t=1, "none", op(t)), t= T)=~C)
    else
        error "polynomial must be linear in the specfied variables"
    fi
end proc
:
p:= randpoly([seq](C[k], k= 0..9), degree= 1, dense);
 p := 74 C[0] + 72 C[1] + 37 C[2] - 23 C[3] + 87 C[4] + 44 C[5]
    + 29 C[6] + 98 C[7] - 23 C[8] + 10 C[9] - 61

CF:= CoeffsOfIndexed(p, C);
 CF := TABLE(sparse, [0 = 74, 1 = 72, 2 = 37, 3 = -23, 4 = 87, 
   5 = 44, 6 = 29, 7 = 98, 9 = 10, 8 = -23, "none" = -61])

CF[0];
                               74

CF["none"];
                              -61

The above works regardless of whether the coefficients or the indices are numeric.

I've had a very similar problem: minus signs appearing as K, plus signs appearing as C, and a few other symbols changed. Just like your case, the problem disappears when the Zoom is greater than 100%. Once it happens to one of my worksheets, all worksheets open in the current session are affected. I haven't had this problem for a few months now, so maybe they're working on it. Anyway, the errors are not saved in the worksheet, nor do they affect the computation; only the display is affected. Closing and reopening my entire Maple session seems to fix it, but closing and reopening individual worksheets does not.

Let's suppose that T and X are two sequences of equal length. For the sake of example, I'll use:

T:= seq(k, k= 1..9):  X:= seq(k^2, k= 1..9):

Then do

TX:= <<T> | <X>>: #make 2-column matrix
plots:-display(
    seq(plot(TX[..k]), k= 1..upperbound(TX)[1]), 
    insequence
);

Then click on the animation controls that appear in the toolbar.

restart:
f:= [
    (x,y)-> -2*x/3+sqrt(x^2+3*y)/3,
    (x,y)-> (-x*y/2+sqrt(x^2*y^2-4*y)/2)*y,
    (x,y)-> 3*sqrt(x*y)
]:
map(
    f-> (
        solve(identity(a*f(a*x,a^p*y)=a^p*f(x,y), a), {p}) 
            assuming a>0
    ), f
);
                  [{p = 2}, {p = -2}, {p = 3}]

subs(x= x-2, convert(subs(x= x+2, 1/(3-x)), FormalPowerSeries));
                        infinity        
                         -----          
                          \             
                           )           k
                          /     (x - 2) 
                         -----          
                         k = 0          

subs(x= x-1, convert(subs(x= x+1, sqrt(x)), FormalPowerSeries));
        infinity                                       
         -----                                         
          \     /      k  (-k)                       k\
           )    |  (-1)  4     factorial(2 k) (x - 1) |
          /     |- -----------------------------------|
         -----  |                   2                 |
         k = 0  \       factorial(k)  (-1 + 2 k)      /

Square brackets cannot be used instead of parentheses. The square of the derivative could be entered as (diff(f(eta),eta))^2 or diff(f(eta),eta)^2; they mean the same thing.

I think that you're causing yourself confusion by making an unnecessary distinction between procedures and subprocedures.  A procedure B only needs to be considered a subprocedure of procedure A if B directly refers to the locals or parameters of A. In the code that follows, SelectX is a subprocedure of ModuleApply, but maxmin and pluralize are not. SelectX needs to be declared inside ModuleApply, but the others do not although they are used within ModuleApply. SelectX itself has a subprocedure, j-> Y[j]=y. Such a procedure is called anonymous because it has no proper name.

The code that you're using as examples uses a very bad programming practice: The procedures return their values through parameters declared evaln rather than through the normal return mechanism.

A group of related procedures and other code can be put together into a module:

optimize:= module()
export
    maxmin:= proc(A::{list, set, table, rtable}(numeric))
    local minv:= infinity, maxv:= -infinity, v;
        for v in A do
            if v < minv then minv:= v fi;
            if v > maxv then maxv:= v fi
        od;
        (maxv, minv)
    end proc,

    Vals:= Record("maxy", "miny", "xmaxs", "xmins")
;
local
    format:= "M%simum y-value is %a and occurs at x-value%s %a.",
    pluralize:= x-> `if`(numelems(x)=1, ["", seq(x)], ["s", x])[],

    ModuleApply:= proc(f, ab::range(realcons), N::posint)
    uses V= Vals;
    local 
        a:= lhs(ab), J:= [$1..N+1],
        X:= a+(rhs(ab)-a)/N*rtable(J-~1, 'datatype'= 'float'), 
        Y:= f~(X),
        SelectX:= y-> [seq](X[select(j-> Y[j]=y, J)])
    ; 
        (V:-xmaxs, V:-xmins):= SelectX~(
            [((V:-maxy, V:-miny):= maxmin(Y))]
        )[];
        plot(
            `[]`~(`[]`~(X,0), `[]`~(X,Y)), 'color'= 'black',
            'thickness'= 0,
            'caption'= sprintf(
                cat(format, "\n", format),
                evalf[4]([
                    "ax", V:-maxy, pluralize(V:-xmaxs),
                    "in", V:-miny, pluralize(V:-xmins)
                ])[]
            ),
            _rest
        )
    end proc
;
end module
:

optimize(x-> 3+10*(-x^2+x^4)*exp(-x^2), -1..4, 150);

eval(optimize:-Vals);
   Record(maxy = 6.08806665291368, miny = 1.39153873739412, 
   xmaxs = [1.63333333333333], 
   xmins = [-0.633333333333333, 0.633333333333333])

optimize:-Vals:-maxy;
                        6.08806665291368

After separating the real and imaginary parts with evalc, as shown by Thomas Richard, there are still a few more things to consider before plotting them. The resulting equations have three (real) variables---p, q, and y---and they can't be explicitly solved for any of those. The command plot is not for plotting equations. An equation in three variables could theoretically be plotted by plots:-implicitplot3d if you specified ranges for all three variables. Given the high degree of your equation, I wouldn't expect very good results from this. If you give a specific value to one variable and specify ranges for the other two, then the resulting equation could be plotted with plots:-implcitplot. This could be done in succession for different specific values of that one variable, producing an animation as in this command:

plots:-animate(plots:-implicitplot, [evalc(a), p= -2..2, q= -2..2], y= 0..2);

The Explore command can also be used in conjunction with plots:-implicitplot, as shown by Kitonum, to produce animations, or they can be used as he showed to produce static plots with a viewer-adjustable free variable (parameter).