The problem of the USB sticks is called an Integer Linear Program or ILP for short. It can be solved like this:

Optimization:-LPSolve(
   32*x + 64*y + 128*z, 
   {15*x + 20*y + 30*z <= 200, x >= 0, y >= 0, z >= 0},
   assume= integer, maximize
);
     [832, [x = 0, y = 1, z = 6]]

The 832 is the maximum number of GB.

There are a vast number of ways to do that. Here's one way:

FullSpeed:= 3600:  #whatever number you want
FullTorque:= 700:  #whatever number you want

Current:= (tq, sp)->
   #Put your equation for current here
   tq^2 * sp #I just made that up
:
P:= plot3d(Current, 0..FullTorque, 0..FullSpeed, grid= [11,11]);
interface(rtablesize= 11): #Allowed Matrix size to print.
Current_Matrix:= op([1,3], P);
#The row index varies the torque; the column, the speed.

Except for one small detail, this is a straightforward elementary textbook linear-programming problem. That detail is that the problem states the resource consumption in units per day, whereas we need days per unit. This is the reason for the 1 /~ in the code below.
 

Download LP.mw

It's a mystery to me why FrequencyTable doesn't offer that option. The similar command Statistics:-TallyInto does offer that option. If you need the output of TallyInto formatted similarly to FrequencyTable, let me know.

The matrix or any other structure or structures that can be assigned to a name can be saved by

save M "MyData.map";

where M is the matrix and the quoted string is any filename of your choice. Then, in another worksheet, do

read "MyData.map":  #Use : to suppress lengthy output!

and then the name M will have the same value as it had in the first worksheet. The above uses a plaintext file. If you read it with a regular text editor, it'll look like an elaborate Maple command starting M:=. If you change the file extension from .map to .m, a faster but non-plaintext file format will be used.

If M contains references to other variables which are not being saved in the same batch, then there can be subtle differences between the .m and plaintext formats. These differences are too deep to explain here.

@Carl Love Here's the code for your most recently described sequence. It uses three auxilliary tables as well as a remember table. One of those tables, used, is declared sparse, which means that the table's entry at unused indices is 0.

restart:
a:= module()
export
   #Record index (position) where values occur, other than most recent.
   pos:= table([0= 0]),
   nextpos:= table([0= 1]), #most-recent position for each value. 
   used:= table(sparse, [0= 2]), #times value has been used
   ModuleApply:= proc(n::nonnegint)
   option remember;
   local p:= thisproc(n-1), r:= `if`(used[p] > 1, pos[p], thisproc(p-1));
      (pos[r], used[r], nextpos[r]):= (nextpos[r], used[r]+1, n); #shuffle down
      r
   end proc
;
   (ModuleApply(0), ModuleApply(1)):= (0,0)
end module
:
seq(a(k), k= 0..30);
0, 0, 0, 1, 0, 2, 0, 4, 1, 3, 0, 6, 2, 5, 0, 10, 3, 9, 1, 8, 4, 
7, 0, 14, 5, 13, 2, 12, 6, 11, 0

Some interesting patterns that I've noticed so far:

1. a(2^n) = 2^(n-2) - 1  (generalized in #12)
2. a(2^n+1) = 3*(2^(n-2)-1) = 3*a(2^n)
3. a(2^n-1) = 3*2^(n-2)-2 = a(2^n+1)+1
4. a(2^n-2) = 0
5. a(2^n+2) = 3*2^(n-4) - 2
6. a(2^n-3) = 3*2^(n-3) - 1
7. a(2^n+3) = 3*2^(n-2) - 4
8. a(2^n-4) = 2^(n-2) - 2 = a(2^n) - 1
9. a(2^n+4) = 2^(n-2) = a(2^n) + 1
10. a(2^n-5) = 3*2^(n-3)
11. a(2^n+5) = 3*2^(n-2) - 5,  n > 3
12. a(m*4^n) = m*4^(n-1) - 1, 
13. a(3*2^n+1) = 2^(n+1) - 3,  n > 1
14. a(3*(2^n-1)) = 2^n - 1
15. a(3*2^n - 1) = 2^(n+1) - 2
16. a(3*2^n - 2) = 0
17. a(4*n+2) = a(n-1)  (generalized in #18)
18. a(n) = a(4^m*(n+2) - 2)

Except where noted (#11 and #13), these apparently hold true for all values of the variables such that the exponents and sequence indices are nonnegative integers. 

While using even a large number of tables is faster than backtracking the sequence, I wonder if anyone can code the above with fewer auxilliary tables.

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.