@_Maxim_ The use of 0 < r < R works in 2D-input only.
In 1D you must write 0 < r, r < R.

Just a comment on the statement by Robert Lopez:
" At one time this quirk could break commands such as dsolve. I don't know if it still does .... ".

Although here is no assignment I tried the following code in Maple 8, Maple 12, Maple 15, and Maple 2017.
It works fine in Maple 15 and Maple 2017.
In Maple 8 neither the symbolic dsolve command nor the numeric worked. In Maple 12 the numeric dsolve worked, but the symbolic didn't.
 

ode:=diff(n(n[1]),n[1])=n(n[1]);
dsolve({ode,n(0)=1});
res:=dsolve({ode,n(0)=1},numeric);
res(1);

The less weird code that follows only failed for the symbolic dsolve command in Maple 8:
 

restart;
sys:=diff(n(t),t)=n[1](t)+n(t),diff(n[1](t),t)=n[1](t)-n(t);
dsolve({sys,n(0)=1,n[1](0)=0});
res:=dsolve({sys,n(0)=1,n[1](0)=0},numeric);
res(1);

The following code containing an assignment to n[2] breaks all of the versions mentioned, numeric as well as symbolic:
 

restart;
n[2]:=47;
ode:=diff(n(t),t)=n(t);
dsolve({ode,n(0)=1});
dsolve({ode,n(0)=1},numeric);

Although a name of a table (here n) only evaluates to the name (n), the command n(t) evalutes to
(table([2 = 47]))(t).

@_Maxim_ This works too:
 

radnormal((sqrt(1-I)*sqrt(1+I)-sqrt(2))*(t+1));

@J F Ogilvie You wrote:
"When I entered this logically correct assignment in Maple 2017.3 (classic interface)

    n := n[1] + n[2] + abs(m) + 1;

.... "
Now what does 'logically correct' mean here? If you are "in the know"  as Robert Lopez calls it then it doesn't make much sense.
He also mentions that this has been with us since Maple V Release 4 (I haven't checked).
Changing it would be another colossal waste of time that could be used elsewhere.
Just stick to using xxx__yyy as vv says.
 

@tsunamiBTP No, the problem lies with constants (like Pi, sqrt(2) or sin(2)) that are not of type numeric. The fraction 1/3 is of type numeric; Pi, sqrt(2), and sin(2) are not.
Actually I mentioned 22/7 as an example for which the substitution Pi=22/7 makes is return true for the original equation. I might as well have used Pi=1/3.
##
Secondly, this "remedy" of replacing Pi with an unassigned name works here, yes, but clearly we could construct examples of equations whose validity depends on Pi being Pi and not just a name.
##
So my point was just to point at the possible cause of the bug.
Basically, I don't think that you should trust is.

 

@tsunamiBTP Your new example doesn't differ in its essence from the simpler one you started with.
The following changes (in order) will give you true:

is(subs(Pi=pi,Y));
##
is(eval(y1 = y2, {Pi = pi, n = 10}));
##
is(eval(y1 = y2, {Pi = pi, n = 10})) assuming real;
##
is(combine(subs(Pi = pi, eq)));
############
You have value in a couple of places. No need for that. You are using sum (not Sum) and sum finds a symbolic sum.
As I was pointing out in my answer, theproblem seems to be the handling of constants that are not rationals like Pi and sqrt(2). I use pi as a name, it might as well have been another unassigned name.

A quite similar question was asked on MaplePrimes here:
https://mapleprimes.com/questions/223417-Physics--Math-With-Maple-DE-Numerically#answer244152

Notice that the OP (Joey Joe) had a dimensionally incorrect version of the Einstein version at first.

@vv I agree with your first 3 lines.
About the last point: I went to the help page
?Atomic variables
and found this gibberish:

Use a literal subscript to form a name that will not be interpreted as an index of any kind. For example, if you want to use both  x  and x[0]  in a worksheet as distinct names, use a literal subscript for
#mrow(msub(mi("a",mathcolor="#AF00AF"),mi("b",mathcolor="#AF00AF"))

Quite simple, right? It made me think that somebody goofed up the help page, but I'm afraid not.

@Muhammad Usman Your solve command in the worksheet HELP.mw actually returns a result.
You have commented out the other attempt to find Sols, thus the graph of P is using that solve command.
But since you have a system with floating point coefficients and the determinant is found to be something like
-3.4266188236700911852*10^(-318)
then it doesn't make sense to convert to rationals.
Consider this simple example:
 

restart;
A:=Matrix([[Pi,1],[sqrt(Pi)*Pi,sqrt(Pi)]]);
LinearAlgebra:-Determinant(A); #0
A2:=evalf(A);
LinearAlgebra:-Determinant(A2); # 0.
A3:=convert(A2,fraction);
LinearAlgebra:-Determinant(A3); # NOT zero

The point is that converting a floating point matrix to rationals often results in a nonsingular matrix although the original really was singular.
##
Instead you can for your system use LeastSquares from the LinearAlgebra package to find the vector x that minimizes A.x-b in the euclidean norm.

LS:=LeastSquares(A,b);
Norm(A.LS-b,2); #Roughly 0.3

and then do
 

res:=eval(theta,var =~convert(LS,list));
P2 := eval(res, t = 0.1e-2);
plot(P2, x = 0 .. R);

@tomleslie Yes, I have the same problem (in Maple 2017.3, Windows 10, 64 bit) that you refer to.
I get

Warning, problem appears to be unbounded
           [0, [xifp = 0, xifq = 0, xp = 0, xq = 0]]

but setting UseHardwareFloats:=false gives me the answer you get.

Since this problem must be solved numerically we need to know the parameters M, R, and alpha, or at least in which ranges they lie.
Even so there is no guarantee of success.

@acer The OP's original code works for me on Windows 10, 64 bit, Maple 2017.3 when using method=interiorpoint.
It doesn't work correctly when using method=activeset.

kernelopts(version);
  Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877

interface(version);
    Standard Worksheet Interface, Maple 2017.3, Windows 10, September 27 2017 Build ID 1265877

@tomleslie You wrote:
" My original post, and all of my subsequent posts correctly solve the problem on 64-bit Windows (as it happen WIndows 7), using LPSolve(). These correct solutions are obtained for Maple 18.02, Maple2015.2, Maple2016.2 and Maple 2017.3. "

Since I understand that the OP (andreasks) works at a university and has to deal with students having different kinds of computers, he must have a solution that works for all types of computers.

As I have posted (and with the default setting of Digits=10 and UseHardwareFloats), your code:

restart;
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       binaryvariables=[x__1, x__2, x__3, x__4, x__5],
                       maximize = true
                      );
#######

does not work on my Windows 10, 64 bit computer in Maple 18.02, 2015.2, 2016.2, and 2017.3.
But it does work (thanks to Markiyan Hirnyk for pointing this out) if Digits is set to 16, thus forcing Maple to use software floats.
Your code as well as the OP's original code also works with default setting of Digits, if we make the assignment UseHardwareFloats:=false: as in
 

restart;
UseHardwareFloats:=false;
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       binaryvariables=[x__1, x__2, x__3, x__4, x__5],
                       maximize = true
                      );
################
restart;
UseHardwareFloats:=false;
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       x__1=0..1, x__2=0..1, x__3=0..1, x__4= 0..1, x__5=0..1,
                       maximize = true
                      );

##############
NOTE 1: I'm not saying that this problem exists on all computers with recent Maple versions (since 18.02) and with Windows 10, 64 bit. Had this been true, my guess is that acer wouldn't have replied as he did in his "specifics" reply.
NOTE 2: The OP's original code works for me with the option method=interiorpoint and with default settings of Digits and UseHardwareFloats (see acer's first reply).

@faisal We need to know what X is, specifically how it is related to f(eta).
Also we must know the values of A and epsilon[1] .
I'm assuming that S = s?

@Markiyan Hirnyk Thanks for the reference.
You wrote in a reply in that link:
 
"I am glad of the workaround by Mariusz. Every big soft includes bugs. In particular, both Maple and Mathematica are not any exception. However, I don't remember a Mathematica's bug which appears in 32-Bit version of OS only."

Actually my system is 64 bit (Windows 10) and raising Digits to 16 works for me:
 

restart;
Digits := 16: 
Optimization:-LPSolve(2*x+5*y, {3*x-y = 1, x-y <= 5}, assume = {integer, nonnegative});

Futhermore, the OP's original code works too:
 

restart;
Digits:=16:
Optimization:-LPSolve(3*x__1+14*x__2+18*x__3+6*x__4+2*x__5, {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10}, x__1 = 0 .. 1, x__2 = 0 .. 1, x__3 = 0 .. 1, x__4 = 0 .. 1, x__5 = 0 .. 1, maximize = true);

Answer:
[29.59999998353600, [x__1 = 0., x__2 = .3999999988240001, x__3 = 1., x__4 = 1., x__5 = 0.]]