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.

It seems that you don't really want plots. What you want is a 50x1000 matrix of computed y values for evenly spaced values of a and x. Each y value can be computed with fsolve, like this

R0:= gamma+2*ln(2)+Re(Psi(1/2+(1+I)*tanh((1+I)*x)/((tau+0.25e-1*a)*y)))+ln(y):
tau:= 9.975:
R0p:= unapply(R0, [a,x]):
R0f:= proc(a,x)
local r:= fsolve(R0p(a,x), y= 0..1);
   `if`(r::float, r, Float(undefined))
end proc:
M:= Matrix(
   (50,1000),
   (i,j)-> R0f(1 + (i-1)*(1-0)/(50-1), 0 + (j-1)*(5-0)/(1000-1)),
   datatype= float[8]
);

This took about a little less than 2 hours on my system, using Digits = 15. It make take less time with a lower value of Digits.

If you actually do want a plot, it will be trivial to generate it from the matrix. It can be a 2D or 3D plot.

eq1:= diff(s(t),t) = (1-phi)*epsilon+(1-rho)*a+(1-f)*alpha*v(t)-(lambda+theta[1]+a+epsilon)*s(t);
eq2:= diff(v(t),t) = phi*epsilon+rho*a+theta[1]*s(t)-((1-f)*alpha+f*theta[2]+a+epsilon)*v(t);
eq3:= diff(e(t),t) = lambda*s(t)-(delta+a+epsilon)*e(t);
eq4:= diff(r(t),t) = eta*i(t)+v(t)*f*theta[2]-(a+epsilon)*r(t);
dsolve({eq1,eq2,eq3,eq4});

The answer is quite lengthy. However, this is a linear system, so the solution is fundamentally simple. If you use dsolve again with numeric values for the parameters, you'll get a simpler solution.

Here is the complete solution done by exhaustive search. My results are identical to what you got with Palisades Evolver, so I'd say that the published paper is wrong! The seven-day problem took 5 seconds; the five-day problem took 41 seconds. How much time did Palisades take?

So people actually make money doing this and they can even publish papers on solution to individual problems? Hmmm, maybe I can find work like that.
 

Download Floral_truck_routing_problem.mw

@Febar You should use = (in various places) instead of :=. Something like this:
 

restart:         Demand not directly related to sales price and second-market demnad Sales returns, consumer's return requests Refurbishing proportion among returned items Profit of manufacturing process     MSR: Integration of manufacturing, selling, and refurbishing processes solve(%, p__N); eval(%, R= r__0 - r__1*b); simplify(%, {d__0 - d__2*b = D__0, r__0 - r__1*b = R});

Download Joint_Pricing_Sol_.mw

This isn't the final simplified form that you were hoping for. But surely you can see that it's impossible to eliminate the variables alpha, g, u, and t using solely the information that you've provided.

Note that my simplify command is in the "simplify with side relations" form that Robert Lopez mentioned.