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!

It is obvious by symmetry^[*1] that the answer to the first two is the same^[*2]. These cannot be done in one step by Maple, but they can be done in two.

f:= (x,y)-> (x+y)*sin(1/x)*sin(1/y):

One-step attempt:

Limit(Limit(f(x,y), x= 0), y= 0); % = value(%);
Error, (in unknown) invalid input: limit expects its 1st argument, expr, to be of type Or(equation, algebraic), but received -abs(y*sin(1/y)) .. abs(y*sin(1/y))

The error message is because the inner limit evaluates to an interval, which indicates oscillatory-yet-bounded divergence. But the outer limit is not prepared to accept an interval as input. We can directly apply the Squeezing Theorem^[*3] like this:

map(limit, limit(f(x,y), x= 0), y= 0);
      0 .. 0

Since the interval has been squeezed to a single point, that point, 0, is the limit.

The syntax for getting Maple to do your third limit is

limit(f(x,y), {x=0, y=0});
 

^[*1]By symmetry, I mean f(x,y) = f(y,x) for all x and y.

^[*2]That the first two limits are necessarily the same makes this exercise a bit weird from a pedagogical perspective: If I were teaching this material, I'd want to highlight that switching the order of the limit operations can make the answers different.

^[*3]Note for experts: MultiSeries:-limit does not allow the possibility of doing this Squeezing Theorem trick because it returns undefined rather than an interval for the inner limit. Perhaps it should be changed to give interval results.

Another way is to use the operator-form of diff, D:

f:= (x,y)-> arctan(x*y^2):
D[1](f)(x,y); D[2](f)(x,y);
 

Your coefficients contain imaginary parts (why?) and thus so does fN(p,t) for almost all real p and t. Those imaginary numbers shouldn't be used in order relations (such as < and >). So, both inequalities in your piecewise FAIL for all your points. Thus the piecewise returns its default otherwise value, 0. It's only a coincidence that 0 was also your cap value; that branch of the piecewise is never being taken.

And why are you using sqrt(t) in your basis functions when t can be negative? Maybe you should try sqrt(abs(t)).