You have a second order ODE. The general solution therefore contains two arbitrary constants. These arbitrary constants can be eliminated if you provide two boundary/initial conditions

Well based on this question (and your previous one), I can't help think that you are in a deep hole and all that I am doing is assisting with the digging!

However if you really, really, really want to do this, then

varnames:=['S','T','U','V','W'];
diff(f(map(entries, varnames, nolist) ), U[0,0]);

ought to work

Rr[x,y] refers to the x-th, y-th entry of the indexable entity Rr.

But it seems that Rr is a (non-indexable) procedure which accepts two arguments!!!!!

So what can the OP mean - once again reduced to guesswork!!!

Perhaps what is required is to apply the procedure Rr() to all entries in an indexable entity - as in the following  'toy' example

  A:=Matrix(2,2,[[1,2],[3,4]]);
  B:=Matrix(2,2,[[5,6],[7,8]]);
  Tr:=proc(x)
                 return x^2:
        end proc;
#
# One way
#
   Rr~(A);
   Rr~(B);
#
# Another way
#
   map(Rr, A);
   map(Rr, B);

On the other hand perhaps what is required is a procedure which takes two values, and can be mapped across two indexable entities as in

  restart;
  A:=Matrix(2,2,[[1,2],[3,4]]);
  B:=Matrix(2,2,[[5,6],[7,8]]);
  Rr:= proc(x,y)
                   return x^2+y^2:
          end proc;
  zip(Rr, A, A);
  zip(Rr, B, B);
  zip(Rr, A, B);

Who knows?

The optimum choice probably depends on what you plan on doing with the columns, once you have them.

  with(LinearAlgebra):
#
# generate a random amtrix, just for illustration
#
  M:=RandomMatrix(10,10);
#
# Generate a list of columns in a few different ways
#
  Column(M, 1..10);
  seq( M[..,j], j=1..10);
  for j from 1 by 1 to 10 do
      M[..,j];
  od;

Don't undertstand the 'for each' construct - but you don't need it when definiing a vector with an embedded function initialiser.

You do however need an explicit value for 'N' in order to define a Vector()

The initialiser function within the Vector definition should be '->' not '-->'

The following works for me

N:=10;
G1:=Vector[column](N,i->f(i,1));
G2:=Vector[column](N,i->f(i,2));

or if you are planning on defining a lot of vectors , then you could use

N:=10;
for j from 1 by 1 to 2 do
    G||j:=Vector[column](N,i->f(i,j));
od;

Since your orginal post stated that you were " simulating the movement of 4 astronomical bodies with Newtons equation", I guessed that this involved solving a system of ODEs numerically. Technically dsolve() returns a procedure (or a list of procedures, or an Array, depending on the options you set in the dsolve() command. Values for all problem variables can be be obtained by appropriate calls to this procedure.

I have added a few execution groups to your worksheet which performs these calls and plots the relevant spacecurves. I have deliberately done this in a step-by-step manner to illustrate the process. I could have achieved the same functionality more efficiently (and with less code to type!), but it would probably be more difficult for you to understand.

Check out the attached

planets.mw

When solve() produces some horrendously long expressions, making it difficult to figure out whihc one(s) might be interesting, it is sometimes useful to use complexplot(), as in

restart;
eqn:= (1/4)*x*(3-2*x)^2/(1-x)^3=18325.73901+exp(36.58420421-10902.09286/T);
plots:-complexplot( [solve(eqn,x)],
                                     T=300..800,
                                      axes=boxed
                                 );

which will plot all the solutions - real and complex. From the output of the above it is farly obvious that over the specified range of the cariable 'T',  one solution lies on the real axis and two don't. If you know Maple's default colour choice for such a plot (darkRed, blue, green) then you can tell immediately that the second solution returned by solve() is the real one. So if the real solution is the one you want then

plot( solve(eqn,x)[2],T=300..800);

will plot it.

You can then delete the complexplot() command from your worksheet to leave

restart;
eqn:= (1/4)*x*(3-2*x)^2/(1-x)^3=18325.73901+exp(36.58420421-10902.09286/T);
plot( solve(eqn,x)[2],T=300..800);

and, hey presto you are a genius!

Well, the fundamental reason you only see one plot is that issue multiple 'restart' commands in your worksheet. Every time you use the 'restart' command you are instructing Maple to clear everything prior to that point!

So the simplest solution would be to delete all the restart commands. However, this means that initialisations/evaluations from the preceding code may influence subsequent code: which happens in your case - and it could take a long time to figure ot precisely why. For this reason (and the fact that your code, as written, is spectacularly inefficient), I have rewritten it so that it generates an Array of plots for each value of the parameter N[c] in your worksheet - this seems to be the only thing which is varied. The only reason I generate an array of plots for each value of N[c] is that your solution process produces multiple real solutions (as well as a lot of complex ones) and I have no on what basis you are selecting a specifi real solution. So for the supplied values of N[c], I generate all real values of 'a', then plot the corresponding polynomial.

Check out the attached - I could probably tidy it up a little more but it is past my bedtime

plotProb.mw

eq1:= x^2+y^2=0.314;
eq2:= y=0.05180967688*x;
ans:=convert(fsolve([eq1,eq2]),list);
map(rhs, ans);

For the two curves cc1, cc2, you need to set the lower limit of the associated plotting range according to fsolve(cc1=P[B]), fsolve(cc2=P[B]). As in

restart;
F1:= (w+u-s)*(-P[A]*P[B]+P[B])+P[A]*(p-s)-c+s:
F2:= P[A]*P[B]*(p-s)-c+s:
c1:= solve(F1, P[A]):
c2:= solve(F2, P[A]):
params := [w=4, u=2, c=2, s=1, p=7]:
cc1:= eval(c1, params):
cc2:= eval(c2, params):
p1:= plot(cc1, P[B]=fsolve(cc1=P[B])..0.5, color=blue):
p2:= plot(cc2, P[B]=fsolve(cc2=P[B])..0.5, color=green):
p3:= plot(P[B], P[B]=0..0.5, color=red):
plots:-display([p1,p2,p3], thickness=2, view=[0..0.5, 0..0.5]);
plots:-display([p1,p2], thickness=2, view=[0..0.5, 0..0.5]);

Tools->Options->Interface

and ensure that the checkbox for "Large toolbox Icons" is ticked

This is probably about as good as you are going to get. Result is acceptable on my 2560X1440 display - but on a 4K, maybe not?

Maple 2015/2016 come with "new" toolbar icons which are supposedly designed to be more readable on "high-resolution" monitors - and these are "more readable" that those for Maple18 on the same 2560X1440 display. But on a 4k display......

This is a regrettable side-effect of the fact that many artefacts (eg icons) in a (eg) a user interface are prgrammed as fixed pixel counts - and as the pisel size reduces, the artefact (aka icon) size reduces. You can make the argument that any applicaion should start with some seriously high resolution icons, then interrogate the display resolution and sample icons appropriately so that they are an accceptable size. I have no argument with this, but I accpet that it would mean a significant rewrite of the UI for many applications.

I feell obliged to point out that this problem is not specific to Maple. I have an issue with several Windows applications, even on a 2560X1440 display .....with a 4k display, ouch :-(

1. Recommended: read the help at ?procedure, and code the required function accordingly
2. Not recommended: use the Matlab-to-Maple translator

I took the code

function R=ff(i,m)
K=1;
Q= zeros( i );
for j=1:2:i
 Q(j)=2*i-K;
 K=4+K;
end
R=Q(m);

and saved it in the file 'test.m'. I then used

restart;
with(Matlab):
FromMFile("D:/Users/TomLeslie/myMaple/test.m", "D:/Users/TomLeslie/myMaple/test.mpl");
read("D:/Users/TomLeslie/myMaple/test.mpl"):

to read the matlab file, output it as a Maple file, then read the Maple file - obviously you will have to change the pathnames in the above to something appropriate for your installation(s). Subsequent invocations of the Maple code, such as

ff(3,1);
ff(19,7);

then returned 5 and 25 respectively. Are these correct? Well maybe!.

I generally find that doing a "machine translation" from one programming language to another is a high-risk activity, so I never recommend it.

However, if you are too lazy to learn how to write simple Maple procedures, then maybe you will be prepared to use the above technique

s is used as a variable in your equations

it is then used as a loop index to provide facilitate your loop statement.

Changing the loop varoable solves your problem, as in

restart:
e1 := 12*g^2+12*h^2+4*i^2+3*j^2=684:
e2 :=  12*l^2+12*m^2+4*n^2+3*o^2=684:
e3 := 12*q^2+12*r^2+4*s^2+3*t^2=172:
e4 := 12*v^2+12*w^2+4*x^2+3*y^2=108:
e5 := 12*g*l+12*h*m+4*n*i+3*j*o=-84:
e6 := 12*g*q+12*h*r+4*s*i+3*j*t=-84:
e7 := 12*g*v+12*h*w+4*x*i+3*j*y=-84:
e8 := 12*l*q+12*m*r+4*n*s+3*o*t=-84:
e9 := 12*l*v+12*m*w+4*n*x+3*o*y=-84:
e10 := 12*q*v+12*r*w+4*s*x+3*y*t=-84:
e11 := g+h+i+j=-1:
e12 := l+m+n+o=-1:
e13 := q+r+s+t=-1:
e14 := v+w+x+y=-1:
e15 := h*i+m*n+3*s*r+4*x*w=-21:
e16 := g*i+l*n+3*s*q+4*x*v=-21:
e17 := i+n+3*s+4*x=-3:
e18 := g*h+m*l+3*q*r+4*w*v=-7:
e19 := h+m+3*r+4*w=-1:
e20 := g+l+3*q+4*v=-1:
e21 := i*j+o*n+3*s*t+4*x*y=-84:
e22 := j*h+o*m+3*r*t+4*y*w=-28:
e23 := j*g+o*l+3*q*t+4*y*v=-28:
e24 := j+o+3*t+4*y=-4:
e25 := j^2+o^2+3*t^2+4*y^2=144:
e26 := i^2+n^2+3*s^2+4*x^2=129:
e27 := h^2+m^2+3*r^2+4*w^2=57:
e28 := g^2+l^2+3*q^2+4*r^2=57:
eqset := {e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15,
          e16, e17, e18, e19, e20, e21, e22, e23, e24, e25, e26, e27, e28};
for zzzz in eqset do print(zzzz) end do;

1. It is not to difficult to write the code to perform what you want. However it looks as if the pdsolve() process is starting to have significant numerical problems.
2. f(x,y) with y=0.02, plotted for x=0..20 starts small, and continues small whilst becoming oscillatory. This characteristic is often associated with numerical noise in the pdsolve() process.
3. the pdsolve() process fails completely for the third PDE (ie when P=10)
4. If pdsolve() can obtain a solution then it is always possible to uses fdiff() to get the functions which you want. But if pdsolve() is generating something close to "numerical noise" under the desired conditions, then performing numerical differentiation on this "numerical noise", whilst possiblel, is not going to be very meaningful
5. Only suggestion I can make at this point is to try to improve the accuracy of the pdsolve process: increase the mesh size, decrease the error tolerances, etc, etc. This is a rather time-consuming process, because solution times will increase a lot, and determinng the optimal options for the pdsolve() command is a bit of an "experimental science".

Anyhow check out the attached, where I have added your current requirements to the end of the worksheet I produced previously. For the reasons stated above, although the code runs, I have little faith in the accuracy of the results being generated

pdeProb2.mw

something like the attached

ruku.mw