The integrand has singularities at -I, 2*I, and 5. All of these are at distance greater than 1 from sqrt(11/2). So no singularities are enclosed by the path of integration. So the answer is 0.

If Maple ever gives you 2D output that you don't understand, just use lprint(%) to see the 1D equivalent. In this case, it will reveal the LerchPhi.

I don't think that dsolve's BVP solvers can handle complex numbers. But your BVP contains a parameter, A3. We can solve for this parameter and then use one of dsolve's IVP solvers, which do handle complex numbers. To solve for A3, note that the right side of your ODE doesn't contain the dependent variable, u. Thus we can use a combination of numeric integration and fsolve.
 

Download Complex_BVP.mw

The cause of the anomaly that you described is that print doesn't return values; it merely displays them in your worksheet. The purpose of print is to display supplementary information from the middle of a procedure. It is not appropriate to use it at the end of a procedure to return the procedure's main return value. If you do use it like this, the actual returned value (not the displayed value) is NULL, which is of type exprseq, which is what you're experiencing. So, all you need to do is change print(Hm) to simply Hm.

I think that you may have the mistaken idea that a "statement" needs to be on a single line. Actually, there's no practical limit to the number of lines that it can take, the number of substatements internal to it, or the depth of their nesting.

Change

if m1 <= 2^(2*k-1) then t1 := (C+(-A1*(m1^2)))/A2 end if;
Do (t1);

to

if m1 <= 2^(2*k-1) then
   t1 := (C+(-A1*(m1^2)))/A2;
   Do (t1)
end if;

And do likewise for the other if-then blocks.

The line

E2:= V(x2);

should be

E2:= V(x2(i));

It's not a source of the error, but I'm concerned about your use of Seed. What is its significance? Why not use rand()?

You should not use randomize() until you have your code running correctly without it because its use will make some bugs harder to find.

Please post code in plaintext form, or post a worksheet. I was unable to test my answer above because I couldn't run your code because it's in a form that can't be copied-and-pasted.

If the function has a Taylor series at x=a, then that's what's returned. For complex-valued functions (which most functions are), having a Taylor series at x=a is equivalent to the function being differentiable throughout an open set containing x=a. Next, if it has a Laurent series at x=a, then that's what's returned. Having a Laurent series at x=a is equivalent to there being a positive integer n such that (x-a)^n*f(x) has a Taylor series at x=a. Finally, if a more-generalized type of series can be found, then that's what's returned. I don't know the details of these more-generalized series, but my first guess is that they're always of the form M((x-a)*g(x-a)), where M(z) is a Maclaurin series and g(x) is a relatively simple function (such as ln(x) or x^r for some r, 0 < |r| < 1) that is not differentiable at x=a.

Here's a somewhat systematic approach to reducing arctrig(algebraic number) to a rational multiple of Pi.
 

Download arctrig_of_algebraic.mw

I consider that last step to be "exact and symbolic" because the min is exactly 0.

The more-typical (but more difficult to implement) solution to this error is to introduce a coefficient called a continuation parameter. The coefficient is usually applied to the more-complicated terms of the ODEs, but it can appear anywhere, and in multiple locations. In this case, I applied it to the nonlinear terms. I used CONT for the coefficient, and changed your first ODE to

ODEforNum1:= 
   diff(u(eta), eta$4) +
   2/(eta+k)*diff(u(eta), eta$3) -
   1/(eta+k)^2*diff(u(eta), eta$2) +
   1/(eta+k)^3*diff(u(eta), eta) +
   CONT*e1*(
      -k/(eta+k)*(diff(u(eta), eta)*diff(u(eta), eta$2)-u(eta)*diff(u(eta), eta$3)) -
      k/(eta+k)^2*(diff(u(eta), eta)^2-u(eta)*diff(u(eta), eta$2)) -
      k/(eta+k)^3*u(eta)*(diff(u(eta), eta))
   ) = 0
;

The only (mathematical) change that I made was the introduction of CONT. The rest is just formatting so that I could easily read the input.

After that, you just need to declare the continuation parameter in the dsolve command:

numsol1:= dsolve(
   {BCSforNum1, BCSforNum2, ODEforNum1, ODEforNum2},
   numeric, output= listprocedure, continuation= CONT
);

The continuation is a kind of bootstrapping. The coefficient starts at 0, which eliminates the complicated terms, and it's gradually increased to 1, which returns the original ODEs. Obviously, in order for this to produce correct results, the parameter can only be used as a coefficient.

To do it with arrow notation:

Vector[row](6, i-> 4*i)

Like this:

plots:-display(
   plot(1/ra, x= 0..0.5, filled),
   plot(1/ra, x= 0.5..1),
   plot([[0.5, 0], [0.5, eval(1/ra, x= 0.5)]]), #vertical line segment
   view= [DEFAULT, 0..10] #Limit the vertical extent.
);

Is that what you mean?

Your subs commands treat y as if it were a function (expression), but it is just a free variable because it is not defined anywhere in your code. Your plot commands treat Psi as if it were a free variable, but it is a function (expression) that you have defined.

Here's a solution that can be easily generalized to higher dimensions and an arbitrary number of bodies.
 

Download Three-body.mw

You misplaced the closing brace. You have

pds:= pdsolve({eq1, eq2, IBC, numeric, spacestep = 1/10});

It should be

pds:= pdsolve({eq1, eq2}, IBC, numeric, spacestep = 1/10);

Then there's a more-subtle problem with your initial/boundary conditions. You have D[1](u)(x,0). Are you sure that that shouldn't be D[2](u)(x,0)? The former can't be handled by pdsolve(..., numeric).

When entering the expression, you must explicity enter the dot:

p:= C.A + A;

This will indicate that you want noncommutative multiplication.