Elementwise operators (those that have ~ appended to their name) behave differently when both of the operands are container objects (such as lists). Observe:

f~([a,b], [c,d]);
       [f(a, c), f(b, d)]

Whereas you were hoping for [f(a, [c,d]), f(b, [c,d])]. To achieve that, you need to use map:

map(f, [a,b], [c,d]);

or, in your case,

map(coeffs, %, [x[1], x[2]]);

The above principle also needs to be applied to the preceeding collect command.

There's possibly an additional complication in this case: coeffs generally returns a sequence as output rather than a list, and you likely need some way to distinguish which coefficients came from which polynomial. In that case, use

map([coeffs], %, [x[1], x[2]]);

Then coefficients of each polynomial will be put in their own list.

Be careful with using the word Vector (with a capital V) because it has special meaning in Maple and may confuse your readers. A Maple Vector is a container object very similar to a list. Your polynomials are in a list, not a Vector. It turns out in this case that had they been in a Vector, the above code would work the same, except that the output would be a Vector.

In 2D input, you need a space between if and the following left parenthesis. This is not needed in 1D input.

Also, the parentheses are not necessary in either form of input.

You can only get a trigonometric solution with real coefficients when the non-real roots of the monic (i.e, M=1) characteristic equation occur as complex-conjugate pairs. That can only happen when C and K (or C/M and K/M) are real.

Perhaps you're interested in something like this:

restart:
f1:= 3.579167 + 3.537500*x1 - 2.645833*x2 - 0.250000*x1*x1 - 0.012500*x2*x1 + 0.175000*x2*x2:
plot(
   [seq(eval(f1, x1= k), k= 3..7)], x2= 3..7, 
   labels= [grape, preference], labeldirections= [horizontal, vertical],
   legend= (ginseng=~ [$3..7])
);

The following at least shows a correct way to input these equations and your known solutions, and plot those solutions using a truncation of the series. I couldn't figure out the pattern of the coefficients for the first solution, so I hard-coded the first 4 coefficients.
 

Download kdv_plot3ds.mw

I don't know if RunWorksheet can do that, and even if it could, it would be a quite convoluted and unusual way to do it.

The easy way is to make your code into a procedure. Then put a call to the procedure into a for-loop, just as you've already suggested. A procedure can have multiple outputs; just separate the outputs with commas.

All function applications in Maple use parentheses, not square brackets. So Re[...] should be changed to Re(...).

You may need to use evalc at the end to achieve the simplification.

If you need further help, you'll need to post a worksheet, not a screenshot.

How are you doing on this problem? I have some hints:

We know from one of the other problems that if f(a) is 1-1 and onto (and thus invertible) then its inverse is f(b) for some b. So the decryption process is essentially the same as the encryption process: Map f(b) over the list of numbers for all possible b. Use convert(..., bytes) to convert the list of numbers into a string T. Then use searchtext("the", T). If this is greater than 0, then T is the decrypted string.

It is difficult, and often impossible, to use an elementwise operator when one of the operands is a container (set, list, table, or rtable) that you want to be treated as a solitary operand.

Here's one more way, which I present as a curiosity rather than as a practical solution:

eval([{1,2}, {2,3}] minus~ S, S= {2})

My favorite way has already been presented by Joe:

map(`minus`, [{1,2}, {2,3}], {2})

Not only is this a bit briefer than

map(u-> u minus {2}, [{1,2}, {2,3}])

it's also more efficient. Keep this is mind whenever you want to map a multi-argument function or binary operator.

You need to define f in your code. From your related Question, "Decrypt Message", I know that f can be defined via

f:= a-> x-> a*x mod 256;

Like this:

for i to 3 do f[i]:= subs(_i= i, x-> P(_i,1,1,x)) od;

Here's a procedure to estimate the order of convergence. The first argument is the converging sequence as a list. The points in the list could be real numbers, complex numbers, vectors, or even more-complicated structures.

OrderOfConvergence:= proc(S::list, {Norm::appliable:= abs}, {diff::appliable:= `-`})
local E:= remove(`=`, map(Norm, map(diff, S[..-2], S[-1])), 0.);
   Statistics:-PowerFit(E[..-2], E[2..])[2]
end proc:

Here's an example of using the procedure to estimate the order of convergence of using Newton's method to approximate sqrt(2):

f:= x-> x - (x^2-2)/(2*x):
Digits:= 200:
S[1]:= 1.:
for k to 9 do S[k+1]:= f(S[k]) od:
ord:= OrderOfConvergence(convert(S,list)):
evalf[10](ord);
                          2.000367876

If the sequence is a sequence of vectors, then the second argument should be LinearAlgebra:-Norm:

OrderOfConvergence(S, LinearAlgebra:-Norm);

0 = ln(1) = ln((-1)*(-1)) = ln(-1) + ln(-1) = 2*I*Pi ?  So, clearly the rule ln(x) + ln(y) = ln(x*y) is not true for all x and all y. However, if the rule is incorrectly applied, the most that you'll be wrong by is an integer multiple of 2*I*Pi. This is mentioned in the Wikipedia article that you linked to, if you read further.

You should show the command that you used for "what I get".

This is pretty close to your desired fieldplot:

plots:-fieldplot([-(x+1)*(x-3)*y, x*(y+2)*(y-1)], x= -1..3, y= -2..1);

Unfortunately, implicitplot3d is one (and I believe the only) plot command from which one cannot extract the data. But your equation can be easily solved for f0 or sigma, each yielding a two-branched solution. That's the way to go.