Current Maple has a form of eval called eval[recurse]. I think that your Maple 15 doesn't have it, but it can be implemented as a one-liner anyway:

EvalRecurse:= proc(e, v::{list, set}(`=`), L:= lhs~(v))
    `if`(has(e,L), thisproc(eval(e,v), v, L), e)
end proc
:

Then, simply do

EvalRecurse(diff(u(Y(xi)), xi$5), {diff(Y(xi), xi) = mu*(1-Y(xi)^2)});

I see three major ways to simplify UV:

1. It's not necessary to assume that variables are real; that is handled automatically by evalc.
2. The number of parameters of most "simple" procedures f can be determined by nops([op](1,eval(f))) . Thus n doesn't need to be passed.
3. By returning a list of two procedures constructed with unapply rather than a list of two expressions you can avoid the need to use possibly contaminated (i.e., already assigned) global variables x[1], y[1], etc. And the results will be easier to use.

Using those, UV becomes

UV:= proc(f::procedure)
local x, y, XY:= map2(index~, [x,y], [$1..nops([op](1,eval(f)))]);
    map(unapply, (evalc@[Re,Im]@f@op@(Complex@op)~)(XY), op~(XY))
end proc
:
#Example usage:
uv:= UV((z1,z2)-> (z1+z2)^2);
  [(x_1, y_1, x_2, y_2)-> -(y_1+y_2)^2+(x_2+x_1)^2, 
   (x_1, y_1, x_2, y_2)-> 2*x_1*y_1+2*x_1*y_2+2*y_1*x_2+2*x_2*y_2
  ]

Complex(uv(1,2,3,4)[]);
    -20 + 48*I

The issue that you're having with the lists Z and F is unrelated to this, but it can also be corrected with unapply, like this:

map(UV@unapply, F, Z);

In this case, it's very easy to express x as an explicit function of y. That would make getting numerical data easier, and the plotting could be done by the more-efficient and more-accurate command plot rather than plots:-implicitplot.

There are several ways to get this explicit function. Given the work already done by @dharr , the easiest way is

F:= unapply(solve(sol2, x), y);

evalf(F(1.1)); #numeric data
plot([F(y), y, y= 0.8..2]);

However, if you were starting from scratch, there are several ways to obtain this explicit function that are easier (from a computation perspective).

1. The FunctionAdvisor 3D plot of ln(z) can be made by

plots:-complexplot3d(
    ln(z), z= -5-5*I..5+5*I, view= [default$2, 0..4], title= ln(z),
    orientation= [-50, 60], labels= [Re, Im, abs@ln](z)
);

Using the real-number 3D plotting command plot3d, the above is equivalent to 

plot3d(
    abs(ln(x+I*y)), x= -5..5, y= -5..5, view= [default$2, 0..4], title= ln(z),
    color= arctan((Im,Re)(ln(x+I*y))),
    orientation= [-50, 60], labels= [Re, Im, abs@ln](z)
);

In complex analysis, for z = x+I*y, the arctan function used in the color option in the plot3d command above is a fundamental function called argument; specifically argument(z) = arctan(y,x). See ?argument, and see ?arctan for the precise definition of arctan(y,x) (essentially it's arctan(y/x) with some small adjustments). Every complex number can be uniquely specified by its absolute value r and argument theta via z = r*exp(I*theta) = r*cos(theta) + I*r*sin(theta). Understanding this equation is one of the most important steps in learning complex analysis.

2. A table of plots can be displayed by plots:-display(A) where A is a Vector or Matrix each of whose entries is a plot (from any plotting command). Examples:

PD:= rcurry(plots:-display, size= [200$2]):
PD(Vector(2, i-> plot(x^i, x= -2..2)));
PD(<seq(plot(x^i, x= -2..2), i= 1..2)>);
PD(Vector[row](4, i-> plot(x^i, x= -2..2)));
PD(<seq(plot(x^i, x= -2..2), i= 1..4)>^+);
PD(`<|>`(seq(plot(x^i, x= -2..2), i= 1..4)));
GT:= GraphTheory:
PD(
    Matrix(
        (4,4), 
        (i,j)-> GT:-DrawGraph(
            GT:-CompleteGraph(i,j), stylesheet= "legacy"
        )
    )
);

1. The 3d plot of ln(z) shown by FunctionAdvisor can be obtained by

plots:-complexplot3d(
    ln(z), z= -5-5*I..5+5*I, view= [default$2, 0..4], 
    orientation= [-45,45], labels= [Re(z), Im(z), abs(ln(z))]
);

2. A table of plots can be displayed via plots:-display(A) where A is a Vector or Matrix each of whose entries is a plot. Examples:

PD:= rcurry(plots:-display, size= [200$2]):
PD(Vector(2, i-> plot(x^i, x= -2..2)));
PD(<seq(plot(x^i, x= -2..2), i= 1..2)>);
PD(Vector[row](4, i-> plot(x^i, x= -2..2)));
PD(<seq(plot(x^i, x= -2..2), i= 1..4)>^+);
PD(`<|>`(seq(plot(x^i, x= -2..2), i= 1..4)));
GT:= GraphTheory:
PD(
    Matrix(
        (4,4), 
        (i,j)-> GT:-DrawGraph(
            GT:-CompleteGraph(i,j), stylesheet= "legacy"
        )
    )
);
 

Regarding your first example, note that even simplify(eq) returns an error. Your integral does not have a name as its 2nd argument; I'm not aware of that being a valid way of specifying an integral. The fact that you've found some simplify calling sequence that does something that you consider useful is I guess just a case of GIGO (garbage in, garbage out).

Looking at the help ?simplify,details, I don't see any possibility of assume= anything being used as the 3rd argument. Nor do I see anything that would allow assume= name::property to be used. The documented usage is only assume= property, which applies the property to all variables in the expression being simplified.

Even in Maple itself, the common use of assign to assign solution values to variables that were previously used symbolically (such as to specify a system of equations as in your case) is a sloppy programming practice that very often causes problems later on. You should use new variables, as in

(x0, y0):= eval([x, y], s)[];

It is not the assign command itself that causes problems; rather, it's assigning values (by any means) to the original variables x and y. Thus, I'd also recommend that Samir's suggestion be changed to

x0:= eval(x, s)  #x0 can be anything other than x or y.

Here is how to use a discrete approximate solution to a numeric BVP. I've used the same BVP as Rouben with a small modification to make it more general by including a symbolic parameter.

restart:
Digits:= 15 #suggested starting value
:
# I've included a symbolic parameter a in the ODE system in order to make this example as general as possible. 
ode:= diff(u(x),x$2) + a*surd(u(x),3)
:
# We need one additional boundary condition for each symbolic parameter:
bcs:= u(0) = 0, u(1) = 0, D(u)(0) = -.01
:
# This procedure solves the system, allowing for an optional discrete approximate solution, 
# and plots the solution.
Solve:= proc({
    DiscAppx::Matrix:= NULL, 
    plotopts::list({name, name=anything}):= 
        ['size'= [300,200], 'numpoints'= 50]
})
local 
    sol:= dsolve(
        {ode,bcs}, 'numeric', 'method'= 'bvp'['midrich'], 
        'maxmesh'= 2^9,
        `if`(DiscAppx=NULL, NULL, 'approxsoln'= DiscAppx),         
        _rest
    )
;
    print(plots:-odeplot(sol, [x, u(x)], plotopts[]));
    sol(0.5) #just to see the format
end proc
:

Solution with no approximation given:

S:= Solve();

S;
  [                                      
  [x = 0.5, u(x) = -0.00278207352306078, 
  [                                      

     d                             -20                       ]
    --- u(x) = -4.47952031890384 10   , a = 0.170377519905245]
     dx                                                      ]

# These names will be the 1st column of the discrete approximation matrix:
lhs~(S);
                     [          d         ]
                     [x, u(x), --- u(x), a]
                     [          dx        ]

npts:= 9: #minimum number of approximation points
X:= Vector[row](npts, i-> (i-1)/(npts-1)):
BuildMatrix:= Model->
    <<lhs~(S)[]> | <X, Model~(X), D(Model)~(X), eval(a,S)~(X)>>
:

Solve(DiscAppx= BuildMatrix(x-> -0.003*sin(2*Pi*x))):

Solve(DiscAppx= BuildMatrix(x-> -0.003*sin(Pi*x))):

Download DiscreteBVPApprox.mw

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?