Your answer for part a) is completely correct.

I agree with Kitonum and VV in this sense: If I were coding these plots for my own purposes, I'd use the evalc method that they show. However, Kitonum and VV seem to think that your trouble is due a lack of Maple knowledge, so they suggest Maple-based solutions. But I can see (with absolute certainty) that you're having trouble understanding the fundamentals of complex arithmetic. My reasons for emphasizing the (x,y) method are pedagogical, to aid your learning.

1. You had already started using it anyway.
2. Your book uses it.
3. It's absolutely essential that you learn how to do it "by hand" for these extremely simple cases.
4. If you were coding a Maple `evalc/...` subprocedure, you'd probably need to use it. You'd at least need to have a fundamental understanding of it.

The reason that problem 3b seems weird is that its author had an additional pedagogical purpose: showing the multiplicative properties of abs and arg:

1. For any complex a and b, abs(a*b) = abs(a)*abs(b). 
2. For complex a and b other than 0, arg(a*b) = arg(a) + arg(b) "mod 2*Pi". In other words, arg has a property very similar to logarithms. This is not a coincidence. (By "mod 2*Pi", I mean that you add or subtract an integer multiple of 2*Pi to keep the result in the correct interval ((-Pi, Pi] in Maple).)

So here's how to apply the (x,y) method to 3b:

arg(2iz) = arg(2i(x+yi)) = arg(2ix + 2i²y) = arg(-2y + 2xi) = arctan(2x, -2y).

abs(4iz) = abs(4i(x+yi)) = abs(-4y + 4xi) = sqrt(16y^2 + 16x^2) = 4abs(z).

I agree with Acer: You should define a scalar operation without using elementwise operators, and then use a single elementwise operation to apply that operation to vectors (or lists).

However, the source of your error is not at all your fault; rather, it seems to be a severe^[*1] bug in the precedence and associativity of the elementwise division operator in 2D Input. A fraction with a long horizontal bar in its 2D representation is meant to have both its numerator and denominator grouped with parentheses in any 1D representation. In other words, this fraction: 

 
                             a + b
                             -----
                             c + d

is a 2D representation of (a+b)/(c+d); this is taught with basic arithmetic in elementary school. But the 2D Input parser is interpreting it as a + b/(c + d).

[*1] The erroneous effect of this bug is severe; fixing it should be trivial.

Yes, go to your Questions page: https://www.mapleprimes.com/users/janhardo/questions
There you will get a list of all your Questions. If you click on the title of any of these Questions, then you can browse it in exactly the same way that you'd browse any other Question on MaplePrimes.

Maple has 5 binary-infix operators, [=, <, <=, .., ::], with the following properties (note that > and >= are always reversed in their internal representation and thus can't be treated as separate, independent operators):

1. They are non-associative.
2. They always have exactly two operands, one on the left and one on the right.
3. These operands can be extracted with lhs and rhs.
4. These operands can be any valid Maple expression, including NULL, a multi-term sequence, or another instance of one of these very same binary-operator expressions. 
5. When the operand is as described in "including ..." in #4, it must be enclosed in parentheses.
6. Expressions formed with these operators are inert.

Thus, equations, inequalities, and combinations thereof with more than two operands can be represented as valid Maple expressions by using parentheses, e.g.,

(1 < x) <= 3
(A = B) = C

The only problem with these is that we don't want those extra parentheses to display in the output. We can achieve this by first using Typesetting:-Typeset on these and then removing the unwanted parentheses, which is why I titled this "Working in the other direction". Here's a procedure that does this:

MultiEquation:= proc(eq)
uses T= Typesetting;
    subsindets(
        T:-Typeset(eq), 
        specfunc(
            And(
                specfunc(T:-mrow), 
                satisfies(
                    f-> membertype(
                        specfunc(
                            identical("=", "<", "≤"), 
                            T:-mo
                        ),
                        f
                    )
                )
            ), 
            T:-mfenced #TypeMK for parentheses
        ),
        op
    )
end proc
:

Examples:
MultiEquation((1 < x) <= 3);
MultiEquation(1 < (x <= 3));
MultiEquation((A=B)=C);
MultiEquation(A=(B=C));
MultiEquation(((A=B)=C)=D);
#Desired parentheses are not removed:
MultiEquation((1 < (x+1)^2) < 2);

As you already realize, each problem can be reduced to a set of equations and/or inequalities in the two real variables x and y. For inequalities, the best command is plots:-inequal.

Here are some mistakes that you made:

a) 1. Note that z + 2 - I = z - (2+I). So, the distance formula is sqrt((x-2)^2 + (y-1)^2).
2. Your sqrt(2) should just be 2.

b) 1. Same mistake as above: distance formula on left side should be sqrt((x-1)^2 + (y-2)^2).
2. Distance formula on right side should be sqrt((x-3)^2 + (y-0)^2).

c) The bar symbol over a complex number or variable is called conjugate in English and Maple. It simply means to negate th​​​​​​e imaginary part. So, the left side is 2*x and the right side is x^2 + y^2.

d) For real x and y, arg(x+I*y) = arctan(y,x), using Maple's two-argument form of arctan.

e) 1. This is an equation not an inequality. Thus, its graph must be a curve (essentially 1-dimensional), not a region (essentially 2-dimnesional).
2. Re(z) = x, Im(z) = y.

f) Same as d.

Opts:= labels= [Re,Im](z), scaling= constrained: #common plot options

(*3a:
For inequalities, use plots:-inequal.
 
You could use the distance formula for abs (absolute value), but it's a
bit easier to use its square.
*)
plots:-inequal(
    [1 <= (x-2)^2 + (y-1)^2, (x-2)^2 + (y-1)^2 <= 4],
    x= 0..4, y= -1..3, Opts
);
(*3b:
This is asking for the points equidistant between two points. Obviously
that's a line. If you solve the abs equation, all the square terms will
cancel leaving just y = x-1. Thus the regular plot command is a good
choice for this one.
*) 
plot(
    solve((x-1)^2 + (y-2)^2 = (x-3)^2 + y^2, y), x= 0..4, y= 0..4, 
    Opts
);
(*3c:
When put in (x,y) form, the equation has both x^2 and y^2. Thus
plots:-implicitplot is a good choice for this one.
*)
plots:-implicitplot(2*x = x^2+y^2, x= 0..2, y= -1..1, Opts);

(*3d: 
For real x and y, argument(x+I*y) = arctan(y,x).
optionsopen is used to highlight the effect of the strict
inequalities.
*)
plots:-inequal(
    [
        1 <= x^2+y^2, x^2+y^2 < 9,
        -Pi/4 <= arctan(y,x), arctan(y,x) < 3*Pi/4
    ],
    x= -3..3, y= -3..3,
    optionsopen= [linestyle= dash, color= red, thickness= 2], Opts
);

Your expression can be displayed in your desired form using this code. This seems elaborate, but the vast majority of this code is to retrieve the mathematical expression (as an ordinary Maple expression) from the displayed form. If the goal was only to get the displayed form without retrieval, it could be done in three lines or so.

restart
:
TS:= module()
uses T= Typesetting;
local
    #Using Typesetting:-_Hold introduces unwanted parentheses. I get
    #around that by using attributes and my own Hold.
    Hold, 
    ModuleApply:= e-> setattribute(T:-Typeset(e), Hold(e))
;
export
    #`macro` is needed because Typesetting:-Typeset declares its 1st
    #parameter uneval.
    Mn:= e-> (parse@sprintf)("macro(%a=%a)", ``||e, T:-mn(""||e)),
        
    Value:= e-> subsindets(
        subsindets(
            subs[eval](``= (()-> args), e), attributed,
            e->
            local p, A:= [attributes](e);
                `if`(membertype(specfunc(Hold), A, 'p'), op(A[p]), e)  
        ),
        specfunc({T:-mn, T:-mi}),
        N-> try parse(op(1,N)) catch: N end try
    )
;
end module
:

Usage:

TS:-Mn(2):
TS(a/(`2`*d)) . ``(b/c);
TS:-Value(%);
lprint(%);

You will need to execute the above to see the output. MaplePrimes will not display it correctly for me.

This Answer is just to provide a clue to the source of the error, not to provide a solution to the problem. The mn is a MathML operator used to make symbols and variables display in an upright style rather than the usual italics. An upright style would ordinarily be used for imaginary i. There are two ways see this in action (in Maple proper rather than MapleFlow, which I don't have):

Typsetting:-mn("i");
                               i

`#mn("i")`;
                               i

The Typesetting commands are only suitable for display, but the second form is just a Maple symbol, so it can be used in computation like any other symbol. For example, you could define

i:= `#mn("i")`:

and then i used in any input expression would display upright in the prettyprinted output.

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