You wrote:

• However, I've been told that maple doesn't run in a distributed-memory parallel sense;

The Grid package implements distributed-memory parallelism.

• as parallelisation in maple is very new.

The Grid package has been around for many years and so has the Threads package (shared-memory parallelism). A related package, process, is so old that it's been deprecated. So it looks like what you've "been told" is very out of date.

I derived a one-line closed-form nonrecursive formula for your sequence:

A_nr:= (n::posint)-> (N-> `if`(n>3*N-2,7,5)*N-n-3)(2^(ilog2(n+1)-1)):
A_nr(0):= 0
:
seq(A_nr(k), k= 0..20);
   0, 1, 2, 4, 3, 6, 5, 10, 9, 8, 7, 14, 13, 12, 11, 22, 21, 20, 19, 18, 17

My previous Answer is still valuable because there are many other sequences that can't be put in a nonrecursive form and require memory of the previously computed values.

I didn't use Maple or any high-powered mathematics to derive the above formula, just basic arithmetic and about 45 minutes of mental perseverance.

It can be done efficiently with two tables, one a remember table (in the background), and one to record the values already used:

a:= module()
local
   a:= table([0= NULL, 1= NULL]), #values already used
   ModuleApply:= proc(n::nonnegint)
   option remember;
   local p:= thisproc(n-1), r:= p-1;
      if assigned(a[r]) then r:= 2*p fi;
      a[r]:= NULL; #Record that r has been used.
      r
   end proc
;
   ModuleApply(0):= 0; #Initialize remember table.
   ModuleApply(1):= 1  #    "         "       "
end module
:
seq(a(k), k=0..20); 
   0, 1, 2, 4, 3, 6, 5, 10, 9, 8, 7, 14, 13, 12, 11, 22, 21, 20, 19, 18, 17

To use NLPSolve, you need to make the objective function into a procedure which numerically evaluates the integral. The minimize command won't help.

f2:= (dat::listlist)-> add(dat[.., 1]) - 2*nops(dat);

or

f2:= proc(dat::listlist) local d; add(d[1] - 2, d= dat) end proc;

or

f2:= (dat::listlist)-> add(dat[..,1] -~ 2);

Note that the syntax dat[.., 1] considers dat as a two-column matrix, the first column being the one that you're interested in.

Each region is convex. The constraints for inclusion in any convex region whose boundary can be described algebraically can be expressed as a simple conjunction of inequalities, using one inequality for each segment of the boundary. This works because any intersection of convex regions is convex. In your case:

R2:= {epsilon <= x_2, x_2 <= y_2, p_B <= y_2, y_2 <= 1}:
R7:= {epsilon <= x_7, x_7 <= 1, y_7 <= x_7, p_l <= y_7, y_7 <= p_B}: 

Whether you use union to express that is immaterial.

I'll assume that you can follow the above to do the rest. If not, let me know.

Note that your diagram presumes epsilon < p_B (strictly). If that's not necessarily true, then these inequalities may require a bit more thought. 
 

Having seen many similar ODE problems, I strongly suspect that you have failed to correctly transcribe from the paper a boundary condition "at infinity". The way that these BVPs are solved numerically in practice is to set a finite value for infinity. For the plot from the paper that you show, I suspect that 6 was used as this "finite infinity", and the condition that you've specifed as (D@@2)(f)(0) = 0 should in fact be (D@@2)(f)(6) = 0. 

Here is the code to produce the plot from the paper. While this code may look complicated, the majority of it is just there to add bells and whistles to the plot and has nothing to do with ODE BVPs:

restart:
ode:= 
   diff(f(eta),eta$3) + diff(f(eta),eta$2)*f(eta) - diff(f(eta),eta)^2 -
   M*diff(f(eta),eta) - A*(diff(f(eta),eta)+eta*diff(f(eta),eta$2)/2) 
   = 0
;
bcs:= 
   f(0)=0, D(f)(0)=1, 
   (D@@2)(f)(inf)=0: #<- This is the critical difference!
params:= A=1, inf=6: #Need to set a finite infinity.
Ms:= [0,1,2,4]:
plots:-display(
   seq(
      plots:-odeplot(
         dsolve(eval({ode, bcs}, {M= m, params}), numeric),
         [eta, diff(f(eta),eta)],
         legend= [M=m],
         color= COLOR(HUE, .85*(m-min(Ms))/(max-min)(Ms))
      ),
      m= Ms
   ),
   size= [800,600], thickness= 2, 
   title= sprintf("\t\t\t\t\t\t\t\t\t A = %3.1f", eval(A, {params})),
   titlefont= [TIMES,BOLD,24],
   labels= [typeset(eta), typeset(`f '`(eta))],
   labelfont= [TIMES,BOLD,16],
   caption= 
      sprintf(
         "\nFig. 2: Velocity profiles for different values of M when A = %3.1f",
         eval(A, {params})
      ),
   xtickmarks= [
      seq(k=k, k= 0..eval(inf, {params})-1), 
      eval(inf, {params})=typeset(infinity)
   ]
);

Yes, this is an easy symbol-for-symbol translation:

Change : to ..
Change = to :=
Change () to [], when the () are used for indexing.

The end result is

dIdQ[nmid-(imax-1)/2..nmid+(imax-1)/2, 1]:= 
   -(omega_t[1..imax, 1] - 
     omega_1[nmid-imax*jmax-(imax-1)/2..nmid-imax*jmax+(imax-1)/2, 1]
    );

I was hoping that someone official from Maplesoft's customer service department would Answer your Question, so I didn't answer earlier. I hope that you haven't given up! I believe that a photocopy of both sides of an official student identification card from an accredited educational institution would be sufficient.

To answer the Commenter, cfdtrader: The price of the software is lower for students.

The units of variance are always the squares of the units of the mean (expected return). So, using your current constraint budget = 10000, the expected return is on the order of a thousand dollars or so, and so the variance is on the order of a million or so. So, your constraint variance <= 0.02 is totally unrealistic and unachievable. You are thinking in terms of squared coeffient of variation (variance/mean^2) for the right side of that constraint whereas the left-side computation is the absolute variance. The former is unitless.

Two solutions:

Use unit budget (b:= 1): If you make the budget = 1, then the solution values for x will represent the fraction of your total budget to invest in each asset, and the solution for the objective function will be the growth rate.

Use coefficient of variation: Do this:

relvar:= 1.; #target relative variance
variance:= expand(X^%T.Q.X - relvar*growth^2) <= 0;

Note that this constraint is equivalent to sigma/abs(mu) <= sqrt(relvar) = relative standard deviation. I think that a target relative standard deviation of sqrt(0.02) = 14% shows a high aversion to risk, but that's a personal decision. Experiment with different values of relvar. In either case, DirectSearch:-Search seems to work better than Optimization:-NLPSolve.

This is a formula that all computer programmers should memorize: For integers a and b, the inclusive number of integers between and including them is abs(b-a) + 1. I call it the fence-post formula because it can be used to answer this: How many fence posts are needed to build a 10-meter fence with posts every 1 meter? Answer: 11. The failure to apply this formula correctly is one of the most common causes of bugs in computer programs. 

In most cases where there is apparently no solution, it is easy to prove that there is indeed no solution.

Proposition: Let p be a prime. Let e[0]=0, and let e[1], ..., e[p-1] be p-1 distinct positive integers such that  p + e[i] is prime for 0 <= i < p. Compute all the remainders e[i] mod p, and suppose that this set of remainders is {0, 1, ..., p-1}. Let q > p be any integer (prime or otherwise). Then the set Q = {q + e[i] | 0 <= i < p} necessarily contains a multiple of p. So, there's a composite in Q.  I urge you as a student of number theory to prove this proposition.

Now let's apply this proposition to the specific cases that you found, 105, 195, 231: The following computations are trivial to do by hand or even mental computation, but let's use Maple:

for n in [105, 195, 231] do
   P:= numtheory:-factorset(n);
   p:= min(P);
   E:= P -~ p;
   R:= irem~(E,p);
   if nops(R) = p then
      printf("Provably no solution for n = %d", n)
   fi
od;
                         P := {3, 5, 7}
                             p := 3
                         E := {0, 2, 4}
                         R := {0, 1, 2}
Provably no solution for n = 105
                        P := {3, 5, 13}
                             p := 3
                        E := {0, 2, 10}
                         R := {0, 1, 2}
Provably no solution for n = 195
                        P := {3, 7, 11}
                             p := 3
                         E := {0, 4, 8}
                         R := {0, 1, 2}
Provably no solution for n = 231

Now, it's probably obvious that for "small" n for which there is no solution, p is most likely to be 3. But it's not always 3. I believe that the first n for which it's not 3 is 95095, for which it's 5. Put 95095 into the for-loop above. Then P = {5, 7, 11, 13, 19}, and for every positive integer e, Q:= P +~ e necessarily contains a composite multiple of 5. Q may also contain a multiple of 3, but for e=6, it does not. 

Given a prime p, the following procedure finds a "no-solution" example whose smallest prime factor is p:

#Find a "no-solution" example with a given smallest prime, p:
FindExample:= proc(p::prime)
local R:= table([0= NULL]), n:= p, p1:= p, r:= 0, k;
   for k to p-1 do
      while assigned(R[r]) do 
         p1:= nextprime(p1);
         r:= irem(p1,p)
      od;
      R[r]:= NULL;
      n:= n*p1
   od;
   n
end proc
:
FindExample~([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]);
[ 6, 105, 95095, 215656441, 5070519693819151, 121621899527020215557, 
  21664759328930363899339259137, 6118163487924060025029213220329793, 
  3042203271558217610793849843552106188681707, 
  410050702193824620789441168165116752238670790985538614407, 
  243521388535217606004348248409375049907078159488916219637351331, 
  65167879497494548612895080562902666093018940639437608943583983945258079733989
  ]

For each of these, there's a good chance that the given example is the smallest such example, but I won't guarantee it.

The values in Q and r are empirical, by which I mean that they come directly from real-world observation (or are made up as a academic example) rather than calculation. Certainly, r is not calculated in the worksheet. The r-values are the expected return rates of the three assets (x, y, z) in the proposed portfolio.

The matrix Q is covariances, not returns (or return rates). It represents that the assets do not grow independently of one another. For example, the negative value (-.20) in the [1,2] position indicates that asset 1 tends to go up when asset 2 goes down and vice versa (covariance matrices are always symmetric).

The expected return (what you're calling EV) is X.r (vector dot product). The author calculates this as add(X[i].r[i], i= 1..3), which is unnecessary because that's exactly the same as X.r.

Your "average" method for EV is invalid (unless perchance the same amount of each asset is purchased, in other words x=y=z).

There is no way to compute r from Q---the concept is not meaningful. It would be analogous to computing the mean of a data set only knowing its standard deviation.

I think that you may have mistaken the letter Q as representing something measured quarterly (i.e., every three months). I don't know why that letter was used, but it has nothing to do with calendar quarters. 

All pictures are rectangles. Your picture is not just a penny; it's a penny within a black square background. That entire square is being squeezed into your plot of a circle. The way around this is to plot an entire square with the values outside the circle being undefined: 

plot3d(
   (x,y)-> `if`(x^2+y^2 <= 1, 0, undefined), -1..1, -1..1, 
   image= "penny.png", glossiness= 0, lightmodel= light4
)

This method is much simpler and much more efficient than those presented in the other two Answers (Kitonum and Tom Leslie):

i:= 1 + ListTools:-BinaryPlace(M[.., j], x, `<=`)

where M is the matrix, j is its sorted column, and x is the target value.