All that you need to do is

collect(vel_qD, diff(s1(t), t))

Since the equations are linear in diff(s1(t), t), you can also use any of the equations that contain diff(s1(t),t) to eliminate it from the other equations. This is easy to do know; so, if you want that, let me know.

Here are your equations in a formatted plaintext form, which is easier to work with (easier for the human; it makes no difference to Maple).

vel_qD := <
   diff(theta1(t), t) = -(diff(s1(t), t))*cos(theta3(t))/(LAD*sin(theta1(t)+theta3(t))), 
   diff(theta2(t), t) = 
     (
        L1*sin(theta1(t)+theta3(t))*diff(s1(t), t) +
        L1*sin(theta1(t)-theta3(t))*diff(s1(t), t) -
        2*LAD*sin(theta1(t)+theta3(t))*diff(s1(t), t)
     ) / 
    (LAD*LBC*cos(theta1(t)+theta3(t)-theta2(t))-LAD*LBC*cos(theta1(t)+theta3(t)+theta2(t))), 
    diff(theta3(t), t) = diff(s1(t), t)*cos(theta1(t))/(L3*sin(theta1(t)+theta3(t))), 
    diff(s2(t), t) = 
       (
           L1*sin(theta1(t)+theta3(t)-theta2(t))*diff(s1(t), t) + 
           L1*sin(-theta3(t)-theta2(t)+theta1(t))*diff(s1(t), t) -
           LAD*sin(theta1(t)+theta3(t)-theta2(t))*diff(s1(t), t) - 
           LAD*sin(theta1(t)+theta3(t)+theta2(t))*diff(s1(t), t)
       ) / (LAD*cos(theta1(t)+theta3(t)-theta2(t))-LAD*cos(theta1(t)+theta3(t)+theta2(t)))
>;

Any member whose reputation is greater than 500 is a moderator. According to https://www.mapleprimes.com/users/mmcdara/reputation, you reached 500 about a month ago, so you've being seeing the red spam icons for about a month, which is why they seem "more and more common" to you.

The red spam icon appears on any unresponded content or any top-level post from any user whose reputation is less than or equal to^[*1] 10. The age of the user's account and their history of prior posting is irrelevant. The purpose of the icon is not to alert the reader that the post is spam; rather, its purpose is to allow you, as a moderator, to decide whether it is spam, which immediately deletes the content and de-activates the user's account^[*2].

As a moderator, you have the power to edit or delete any content whatsoever^[*3]; of course, you're honor-bound to use this power with "moderation" and discretion^[*4]. But there's a distinction between simply deleting a post and deleting it as spam---the latter deletion de-activates the user's account as well.

Spam, as a noun in it's original Internet-related usage, refers to content that is totally irrelevant to the forum where it's posted, not content that is merely and unintentionally irrelevant to the thread in which it's posted. Usually, it's advertising. On MaplePrimes, I've deleted everything from ads for plumbers in Toronto to whores in Mumbai. There's also a newer usage of spam as a verb to refer to the frequent posting of low-quality annoying content (which is nonetheless still marginally relevant content). I'm not a fan of this new usage, but unfortunately it seems well established. The moral gulf between the two activities is too wide to risk conflating them ("to tar them with the same brush" as the English saying goes): The former is highly unethical and always intentional, whereas the latter is merely rude and often unintentionally so.
 

^[*1]The MaplePrimes documentation implies that the criterion is (reputation < 10). This is incorrect! In reality, it's (reputation <= 10). I suspect that this was a programming error and that the lower criterion was the intention.

^[*2]We're discovered over the years that it's essential that a spammer's account be de-activated. If this is not done, they post it faster than it can be deleted. In the rare case that something/someone gets deleted inappropriately, the adminstrators can restore it. So, don't worry about that. The important thing is that you delete spam as soon as you see it. This site used to be inundated with spam back in the days when moderators could merely delete it.

^[*3]If content has been responded to, it may be necessary to first delete the responses.

^[*4]I believe that moderators should keep in mind at least a vague notion of some ethically based principles of editorship such as we all should have learned in our formal educations.

The computations (of the level curves) for a contour plot are essentially three-dimensional, even if the final result happens to be displayed as a two-dimensional projection, so a two-dimensional coordinate transformation (such as polar) does not apply. You must use a three-dimensional transformation. The rest of what I would say has already been said by Thomas Richard.

Regardless of whether the range is numeric, a command to do it is

subsindets(eqn, specfunc({sum,Sum}), S-> op(1,S))

where eqn is the equation that you refer to. This works because the summand is operand 1 of each sum or Sum. 

You should never copy-and-paste the output of one command to the input of another. That way madness lies. 

The command codegen:-GRADIENT does what you want:

codegen:-GRADIENT((x,y)-> (x^2-y^2));

See also codegen:-JACOBIAN and codegen:-HESSIAN.

General advice:

Don't use the with command if you don't need to. Instead, use package prefixes (A:-B is command B from package A). Almost everything that can be done with with can be done this way instead. In particular, with(VectorCalculus), literally changes the definitions of numerous (about a dozen) fundamental Maple commands, including D, diff, and even + and *.

Don't use assume if you don't need to. Instead, use assuming.

Since you're using an integral transform, and integration doesn't care about removable discontinuities, why use piecewise? Doesn't this do what you want?

with(inttrans):
sinc:= x-> sin(Pi*x)/x/Pi:
SE := x->sinc(x)*exp(-x^2):
plot(abs(fourier(SE(x),x,s)),s=0..5)

One problem with your most-recent variation is that you've made the coefficients of the integrals 5/Pi. The coefficient does not depend on the number of term. It's always 2/P, where is the period of the overall function.

Here's a little procedure that returns the coefficients, the sum, and/or the plot, depending on the output options that you choose.

FourierSum:= proc(
   f, Per::range, n::nonnegint, 
   {output::{identical(sum,coeffs,plot), list(identical(sum,coeffs,plot))}:= sum}
)
local 
   k, P:= (rhs-lhs)(Per), CS:= [cos,sin], cs, x, X:= 2*Pi*x/P,
   Out:= table(), Output:= `if`(output::list, x-> x, `[]`)(output)
; 
   Out[coeffs]:= [seq(<seq(2/P*int(f(x)*cs(X*k), x= Per), cs= CS)>, k= 0..n)];
   if ormap(`in`, {sum, plot}, Output) then
      Out[sum]:= unapply(add(Out[coeffs] .~ [<1/2,0>, seq(<CS(X*k)[]>, k= 1..n)]), x);
      if plot in Output then 
         Out[plot]:= plot([f,Out[sum]], Per, _rest);
         if nops(Output) > 1 then userinfo(1, FourierSum, NoName, print(Out[plot])) fi
      fi
   fi;
   eval(`if`(nops(Output)=1, op, `<,>`@op)(Output), op(eval(Out)))  
end proc:

Usage:
 
infolevel[FourierSum]:= 1:
f:= t-> t^2:
FourierSum(f, 0..2*Pi, 5, output= [coeffs, sum, plot]);

Your D commands are correct. Why do you doubt it? What output do you get?

While I would just use @@ like the other responders, I thought that you should see a better way to do it as a simple loop (in the hopes that this will improve your coding):

App:= proc(F, n::nonnegint)
local r:= _rest;
   to n do r:= F(r) od;
   r
end proc:

Notes:

The keyword _rest represents arguments that are passed but not declared. This allows the code to iterate procedures F of any number of arguments, including 0.

Only one temporary variable is needed; no swapping (your z1:= z2) is needed.

If the index variable (your j) isn't used in the loop, then the for j from 1 can be omitted.

As noted by the others, the zeroth iterate of a function is the identity function, not the null function.

The problem of redundant calls to the procedure that evaluates the derivative can be solved by something called a remember table. This detects when the procedure is called with arguments that have previously been used and immediately returns the remembered result. To get a remember table in this case, you just need to add the option proc_options= remember to your unapply call. I implemented this in the code below. I also removed numerous other redundancies.

restart:

a:= 0; b:= 1; eps:= 1e-3:
f:= unapply(2*x*(x^2+y), [x,y], proc_options= remember);
G:= simplify(dsolve({diff(y(x),x)=f(x,y(x)),y(a)=1}));                    
N:= 15: h:= evalf((b-a)/N):
X:= Array(0..N, i-> a+i*h, datatype= float):
Y:= Array(0..N, datatype= float);
eps:= 29*eps;
Y[0]:= 1;
(h2, h3, h43, h6):= (h*[1/2, 1/3, 4/3, 1/6])[]: 
for i from 0 to 2 do
   X_i:= X[i]:  Y_i:= Y[i]: 
   t1:= evalf(f(X_i,    Y_i)):
   t2:= evalf(f(X_i+h2, Y_i+t1*h2)): 
   t3:= evalf(f(X_i+h2, Y_i+t2*h2)):
   t4:= evalf(f(X_i+h,  Y_i+t3*h)):
   Y[i+1]:= evalf(Y[i]+h6*(t1+2*t2+2*t3+t4)):
end do;
(S_im1, S_i):= (Y_i, Y[3]):
for i from 3 to N-1 do 
   Y[i+1]:= evalf(Y[i-3]+h43*(2*f(X[i],Y[i])-f(X[i-1],Y[i-1])+2*f(X[i-2],Y[i-2]))):
   S_ip1:= evalf(S_im1+h3*(f(X[i+1],Y[i+1])+4*f(X[i],S_i)+f(X[i-1],S_im1))):
   if abs(Y[i+1]-S_ip1) >= eps then Y[i]:= S_i end if:
   (S_im1, S_i):= (S_i, S_ip1): 
end do;

plot([rhs(G), <X^+|Y^+>], x= a..b, color= [yellow,black], style= [line,point]);

The number of entries in a procedure's remember table is the number of times that the procedure's code was executed. This can be counted by

numelems(op(4, eval(f)));

In this case, it's 37. One of those is the call to the symbolic solver.

If you look closely at Drumeq||1, for example, by simply entering

Drumeq1;

at a command prompt, you'll see that it is a sum of both scalar and Matrix terms. I doubt that this was your intention, and even if it was, I seriously doubt that CodeGeneration is prepared to handle such a highly symbolic concept. (The concept could possibly be meaningful if the intention is that symbols in the scalar part will eventually have Matrices substituted for them.)

How about

evalc(exp(I*2*%Pi/3));

which returns the intermediate form that you want, followed by 

value(%);

to get the final form?

Any name can be made inert by preceding it with % (including names that are already inert or begin with %), and each call to value removes one level of inertness.

In Maple, invocations of named functions, such as sin, always require parentheses. Maple will interpret siny as just another unassigned name.

Once you've entered the ODE correctly, you'll need to specify numeric to get any solution at all from dsolve, because it has no symbolic solution for this one.

Once you have a solution from dsolve(..., numeric), it's possible to plot it with plot, but easier with plots:-odeplot.

restart:
de:= diff(y(x),x) + y(x) = 2*x*sin(y(x)):
cond:= y(0)=1:
y1:= dsolve({de, cond}, y(x), 'numeric'):
plots:-odeplot(y1, [x, y(x)], 0..10);

To disambiguate, use ':-area'. The :- prefix (with no left operand) means "Use the global version of the name rather than the version from any module or package."

I wholeheartedly agree with Robert Lopez about the status of the student package. However, its use is not the cause of your error. The cause is the order in which you've specified the variables of integration. They should go from the innermost integral on the left to the outermost integral on the right. In other words, place the x= ..., y= ..., z= ... in the same order as you would place the dx, dy, dz if you were writing the integral by hand. Note carefully that this is not the same order that you would write the limits themselves if you were writing them on the integral signs!