Here is a way completely different from Kitonum's. In the procedure below, I take the two quadratic solutions returned by solve and determine which is the positive branch numeric evaluation.

PositiveRoot:= proc(
     eqn::{algebraic,`=`},
     var::name,
     {evals::set(name= complexcons):= indets(eqn, And(name, Not(constant))) =~ 0}
)
local
     Sol:= [solve(eqn, var)],
     E:= evalf(eval(Sols, evals))
;
     if nops(E) = 0 then [][]
     elif nops(E) = 1 then Sol[]
     elif nops(E) > 2 then
          error "Procedure only intended for quadratic equations."
     elif Im(E[1]) > Im(E[2]) then Sol[1]
     elif Im(E[1]) = 0 then
          `if`(E[1] > E[2], Sol[1], Sol[2])
     else
          Sol[2]
     end if
end proc:

Example of use:

PositiveRoot(eq1, R);

I echo Markiyan's concerns: The extremely large coefficients, especially in eq3 and eq4, overwhelm the significance of the relatively small final results. Consider using exact arithmetic.

The proximate cause of your error message is that you enclosed the boundary conditions in square brackets. Thus, you have two sets of brackets in the call to pdsolve. 

Any finite, simple, closed region of a computer plot is a polygon, being composed of a finite number of line segments, however small. Thus it can be filled with plots:-polygonplot. You may need to extract the vertices from the original plot structure to do this. If you post your code, I can show you how. 

Indeed A^%T is not even necessarily square, let alone invertible. Use command LinearAlgebra:-LeastSquares(A, b, optimize). You can also multiply b by the pseudoinverse of A, which you can obtain as LinearAlgebra:-MatrixInverse(A, method= pseudo).

Your syntax is perfect, and it gives me the answer -3/2. Try using the restart command.

Taking the constants to be 1 and plotting shows that there is only one positive solution. The solve gives me one answer plus the warning. So, just ignore the warning and use the answer, which is h*c/T/k/(LambertW(-5*exp(-5) )+5).

I'm not familiar with the coeffs option to dsolve. Where is it documented? If I try it, I just get an error message that leads me to believe that no such option exists. But just trying without any options, I get an answer in terms of integrals of Airy functions. 

Apply the command fnormal to the results of floating-point computation before making an equality comparison with zero. 

plots:-display([
     plot(p),
     plots:-textplot(
          zip((a,b)-> [a[], ""||b], p, [$1..nops(p)]),
          align= {above, left}
     )
]);

If you cancel theta[m]s in your x, you get

From this, it is clear that x > 1/4 for the range of phi that you are using. So that plot showing x going to 0 is bogus. (Note: x is the vertical axis---labelled "Cooperation level.")

Where are you reading the "description of the function" that you mention? Is it at ?IsSimilar? What version of Maple? Looking in Maple 16 help, I see no mention of a third argument P. My help says that the third argument can be output= 'C' (the C being literal). Thus, to get P, you can do

P:= LinearAlgebra:-IsSimilar(A, B, output= 'C');

No initialization of P is required.

nops does not work on tables. Use numelems instead.

This procedure is not exactly elegant, but it seems a difficult problem, and this procedure solves it in full generality. It will handle any expression with products that have any number of factors that are Sums. It generalizes from two to any number by applying subsindets recursively. It generates a unique global index name for every distinct Sum. The sequence of indices defaults to _j1, _j2, ..., and the default prefix _j change be changed by using the indexprefix option to the procedure.

ProductOfSumsToSumOfProducts:= proc(expr, {indexprefix::name:= :-_j})
global _indexsuffix;
local
     r, #return value

     #extracts the factors of a product that are Sums
     Extractor:= P-> select(type, {op(P)}, specfunc(anything, Sum)),

     #selects products having multiple Sum factors
     Selector:= And(`*`, satisfies(P-> nops(Extractor(P)) > 1)),

     #transforms a product containing two or more Sums by combining
     #an arbitrary two of them
     Transformer:= proc(P::`*`)
     local
          S:= Extractor(P), S1, S2, S12,
          j1:= indexprefix||(:-_indexsuffix+1), #1st index
          j2:= indexprefix||(:-_indexsuffix+2)  #2nd index
     ;
          :-_indexsuffix:= :-_indexsuffix+2;
          S1,S2:= S[1..2][]; #arbitrary 1st two
          #combine them
          S12:= Sum(
               Sum(
                    subs(op([2,1],S1)= j1, op(1,S1)) * subs(op([2,1],S2)= j2, op(1,S2)),
                    j1= op([2,2],S1)
               ), j2= op([2,2],S2)
          );
          algsubs(S1*S2= S12, P)
     end proc
;
     if assigned(:-_indexsuffix) then  unprotect(':-_indexsuffix')  else  :-_indexsuffix:= 0  end if;
     r:= expr;
     while indets(r, Selector) <> {} do  r:= subsindets(r, Selector, Transformer)  end do;
     protect(':-_indexsuffix');
     r
end proc:

Example of use:

e:= Sum(x[k],k=1..n)*Sum(y[k],k=1..m)*Sum(z[k],k=1..p);

ProductOfSumsToSumOfProducts(e, indexprefix= j);

Here's a standard technique that works in almost any language: Put your program in a while loop. After your program, but still in the loop, get the user's input regarding continuing. Like this:

UserProgram:= proc()
local Again:= Y;
     while Again in {Y,y} do

          # Your original program here
          # ...
          # end of your original program
       
          Again:= "";
          while not Again in {Y,y,N,n} do
               Again:= readstat(cat(Again, "Run again? (Y/N) "));
               if not Again in {Y,y,N,n} then  Again:= "Invalid input. " end if
          end do
     end do
end proc;

That being said, I find the dialog boxes produced by readstat (and the whole family of a commands for taking user input) very awkward. If you want to write Maple programs that process user input from the keyboard, use a Maplet. But that's like learning a whole new language.

If your user wants to run the program again, they can move the mouse pointer up to the program call, change the parameters, and press Enter.

Here are two ways that you can use plots:-display for this. The first prints the plot one after the other. The second uses the fancy Array form of display, and a display inside the display, to display the plots horizontally. The only library command that can do this is display.

Download Array_plot.mw