Do you want to find one feasible point for every possible setting of your variables? Then you need to call LPSolve multiple times anyway. If you have 11 of these "boolean" (two-valued) variables, that's only 2^11 = 2048 calls to LPSolve. In my opinion, that's trivial compared to the huge performance and accuracy penalties you'd pay for making the constraints nonlinear by using a1, a2, etc., as decision variables.

Note that there's no need to use a nested for loop to set multiple variables to multiple values. There are numerous ways to iterate over a Cartesian product that have been discussed many times on MaplePrimes. For example, your double for loop above could be replaced with

nv:= 2: #Number of "boolean" variables.
Cons:= {x1 >= 0, x2 >= 0, x1 <= 1000, x2 <= 1000, a1*x1+x2 <= 700, x1+a2*x2 >= 300}:
BV:= combinat:-cartprod([[0,1]$nv]): #or [[-1,1]$nv]
while not BV[finished] do
     A:= BV[nextvalue]();
     Sol[A[]]:= Optimization:-LPSolve(0, eval(Cons, [a||(1..nv)]=~ A))
end do:

op(eval(Sol));

There are a few errors in your procedure. The first is my fault, and also affects quadsum. I erroneously thought that isqrt(n) returned essentially floor(sqrt(n)) and that iroot(n,r) returned essentially floor(n^(1/r)). They both return something closer to round(...) rather than floor(...). So, replace iroot with this:

Iroot:= (n::nonnegint, r::posint)-> trunc(evalf(n^(1/r))):

That should be good enough for the smallish integers that we're working with here.

The next error is that iquo(n,3) needs to be iquo(n,2). If there's a solution of the form x=y, then 2*x^3 = n, not 3*x^3. Likewise, the check 3*x3 <> n needs to be 2*x3 <> n.

Finally, x=x+1 needs to be changed to x:= x+1.

It is not an error, but there's no need for listcub to be global. I made it local.

So, the fully corrected procedure (with Iroot above) is

cubesum:= proc(n::nonnegint)
local
     k:= 0, listcub:= table(),
     x:= Iroot(iquo(n,2),3), y:= x, x3:= x^3, y3:= x3
;
     if 2*x3 <> n then  x:= x+1;  x3:= x^3;  y:= x;  y3:= x3  end if;
     while x3 <= n do
          y:= Iroot(n-x3, 3);  y3:= y^3;
          if x3+y3 = n then  k:= k+1;  listcub[k]:= [y,x] end if;
          x:= x+1;  x3:= x^3
     end do;
     convert(listcub,list)
end proc:

cubesum(1729);

Just use and between the two conditions and drop the second if.

if irem(cubeprod,3)=0  and modp(cubeprod,2)<>0 then
     cubesum(cubeprod)
end if;

You claim that your second form, the nested if statements, is also erroneous; however, I can find no error in it. What error does it give you? It will do the same thing as the the form with and. The and form is easier to read. 

I notice that you use both irem and modp. This is confusing because they do the same thing. For pure integer arithmethic, use irem; it is a little bit faster because modp also processes matrices and polynomials.

You don't need "to plot inside the for loop." You need to plot after the for loop is finished, like this:

plot(zip(`[]`, vds, convert(Qg, list)));

I don't understand your second question "How can I set the x-axis of the plot in decimal notation?" The above command will give you an x-axis from 0 to 2. Isn't that what you want?

And please, please try to present your code in a neater way. Some tips:

1. You can upload an actual worksheet.
2. The space bar is your friend.
3. Indent continuation lines.
4. Make sure that input and output are clearly distinguished. For example, I always use bold for input.
5. Don't include > at the beginning of your input lines. People need to manually delete those before they can execute your code.
6. Most output is not worth including in a post.
7. Select format "preformatted" from the MaplePrimes toolbar before pasting monospaced output. If you don't, the exponents all get pushed to the left, as you can see above.

Regardless of differences between versions, I can make the whole thing run in 12 seconds by a trivial application of evalhf. After your definition of procedure J, add the line

Jhf:= (x,y)-> evalhf(J(x,y)):

Then plot3d(Jhf, ...) instead of plot3d(J, ...).

I thought that plot3d already used evalhf by default, but I guess that I was wrong about that.

One benefit of the OP having posted a screen shot rather than plain text is that it clearly shows that both P and P appear in the input. My guess is that they are being treated as different variables and that this is the source of the anomaly. Furthermore, I guess that the upright P is actually the Greek letter capital Rho copied from a palette.

Simply include the argument scene= [x(t),y(t)] in your DEplot command and you'll get this:

Download RiemannSum.mw

While Kitonum's brute-force search does find the solutions to the equation, it doesn't use the algorithm given in the Question, which is a much better algorithm. I implemented the algorithm in the Question, and then I made a few improvements to it.

My procedure, like the algorithm in the Question, does not list the symmetric and negative solutions. If those are desired, they are trivial to generate by post-processing the output of the procedure.

quadsum:= proc(n::nonnegint)
local
     k:= 0, mylist:= table(),
     x:= isqrt(iquo(n,2)), y:= x, x2:= x^2, y2:= y^2
;
     if 2*x2 <> n then   x:= x+1;  x2:= x2+2*x-1;  y:= x;  y2:= x2  end if;
     while x2 <= n do
          y:= isqrt(n-x2);  y2:= y^2;
          if x2+y2 = n then  k:= k+1;  mylist[k]:= [x,y]  end if;
          x:= x+1;  x2:= x2+2*x-1
     end do;
     convert(mylist, list)
end proc:

Applying this procedure to random input on the order of 100 billion (10^11), it returns the solutions nearly instantaneously (156 ms). I let Kitonum's procedure run for over 23 minutes with this input with no results.

n:= eval('rand'(2^17..2^18)()^2 + rand(2^17..2^18)()^2);

st:= time():  quadsum(n);  T:= time()-st;

[[221539, 205782], [236579, 188298], [238766, 185517],
  [257901, 157838], [260634, 153283], [269181, 137722],
  [270774, 134563], [281011, 111618], [285826, 98637],
  [287843, 92586], [297357, 54814], [297549, 53762],
  [297978, 51331], [300477, 33754], [302242, 8691], [302323, 5154]
  ]
                             0.156

The algorithm used by isolve is even faster, by about a factor of two.

Sin:= unapply(subs(infinity= n-1, convert(sin(x), FormalPowerSeries)), [x,n]);

plot(Sin(x,20), x= -2*Pi..2*Pi);

Yes, plots:-arrow (not plottools:-arrow) is the answer. Here's a procedure to do what you want. Any optional arguments that you pass to this procedure (such as options controlling the size, shape, and color of the arrows or general plot options such as those controlling the display of axes) will be passed on to plots:-arrow.

A:= [<1,1>, <2,2>, <3,3>]:

HeadToTail:= proc(L::list(Vector(2)))
local L0:= [<0,0>, L[..-2][]];
     plots:-arrow(L0, L -~ L0, args[2..]);
end proc:
     
HeadToTail(A);

I entered your equation into dsolve in Maple 16 and I got a solution.

eq:= (-2*r*beta(t)+J)*diff(beta(t),t$2)+sin(beta(t))*(r*diff(beta(t),t)^2+g)=0:
dsolve(eq);

In this case, it means that you made a syntax error. The y in your difeq should be y(x). That is how you indicate to Maple that y is the dependent variable and x is the independent variable.

Of course you get errors: You're using Matlab code in Maple. The syntax is completely different. Maple syntax is much more compact and intuitive. Here is how to do it in Maple:

a0:= 80;  a1:= -5;  b1:= -5*(3)^0.5;  w:= Pi/12;  k:= 0.2;
U0:= [65,70,95]:
T:= 0..100:
c0:= u0-a0-(k^2*a1-k*w*b1)/(k^2+w^2);
c1:= (k^2*a1-k*w*b1)/(k^2+w^2);
d1:= (k*w*a1+k^2*b1)/(k^2+w^2);
u:= a0+c0*exp(-k*t)+c1*cos(w*t)+d1*sin(w*t);
plot(
     [seq(u, u0= U0)], t= T, legend= [seq('u'(0)= u0, u0= U0)],
     color= [red, black, green], linestyle= [solid, dash, solid]
);

Using the S and t from your original Question, the idea of Tom Leslie can be implemented like this:

S3:= zip((x,y)-> [x[1],x[2]*y], S, CurveFitting:-ArrayInterpolation(rtable(t), S[..,1]));

By the way, "serie" is not an English word. The word is "series" whether it's singular or plural.