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)).

One way to answer the question is to solve the linear system. Since the system is homogenous (right sides are all 0), this  can be done like this:

A:= <
   1,  1,  1;
   1, -2, -1;
   2,  1,  3;
  -1,  2,  1;
>:
nops(LinearAlgebra:-NullSpace(A));
                               0

The nops is to count the vectors in the basis of the null space, which is the same thing as the dimension of the solution space.

I don't understand Kitonum's Answer at all, except that I agree with his first sentence "The dimension of the solution of a system is equal to the number of arbitrary constants in its solution." In particular, I don't see what t and derivatives have to do with this problem.

The problem can be solved with just fundamental arithmetic commands. No "solve" type command is needed. And, egads, certainly no nested loop is needed either.

Sol:= b-> `if`(igcd(b,10000)=1, 2391/b mod 10000, FAIL):

Your example:

Sol(297);
      9503

You could easily use the above to generate a list of all solution pairs, but I don't see the point of doing that. A solution exists iff b and 10000 are relatively prime, in which case the solution is unique (modulo 10000).

To answer your titular question, "Can isolve be used with irem?": The msolve command would almost always be a better choice than doing that. But for the case at hand, even using msolve is overkill, and inefficient.
 

I find the syntax of implicitdiff overwhelming, and I always need to try several things before I find what works.

I will introduce a new dependent variable, F, to represent the expression that you want to differentiate. I'll consider c and r also as dependent variables. I'll consider the independent variables to be beta, varphi, and x, in that order. (I don't think that it actually matters whether we consider x to be a variable or a constant.) So, the input system and variable specifications are

restart:

Sys:= {
   F = 
      c*(r-1)*exp(x*beta)/((1+varphi*exp(x*beta))
      *(-varphi*exp(x*beta)*r+varphi*exp(x*beta)+1)),
   c = exp(exp(x*beta)*(r-1)/(1+varphi*exp(x*beta))),
   ln(r) = varphi*exp(x*beta)*(r-1)/(1+varphi*exp(x*beta))-1
};
Vars:= {F, c, r}(beta, varphi, x);

It's not clear whether you want the first derivatives with respect to beta and varphi or the mixed second derivative with respect to both, or all three of those. So, I'll do all three:

implicitdiff(Sys, Vars, {F}, beta);
implicitdiff(Sys, Vars, {F}, varphi);
implicitdiff(Sys, Vars, {F}, beta, varphi);

Since I hate redundant code, here's a way to combine those three commands (although I admit that this requires expert-level Maple knowledge):

(op@curry(implicitdiff, Sys, Vars, {F})@op)~([[beta], [varphi], [beta,varphi]]);

It's a great Question about a not-well-known feature of Maple. The syntax required to use a procedure's indices as if they were arguments is awkward, but still doable. The following will work if you're always going to call f with an index or indices:

f:= (x,y)-> a[op(procname)]*x + b[op(procname)]*y;

Note that there's no index on f in the definition, so it's impossible to use indices as if they were named parameters, such as your i; you must use the keyword procname.

The proper call to f is just as you suggested, f[3](x,y). But if you now call f without indices, you'll get a nonsense response, although not an error message. So, here's a more-robust treatment, essentially making the index a type-checked parameter:

f:= proc(x,y)
local i;
   if procname::'indexed' then 
      i:= [op(procname)];
      if nops(i) <> 1 then error "Exactly 1 index required" fi;
      i:= i[];
      if not i::'nonnegint' then error "Index must be nonnegative integer" fi
   else
      error "An index is required"
   fi;
  
   a[i]*x + b[i]*y
end proc:

Of course, it's up to you to decide what types and numbers of indices are allowed. Maple imposes no restrictions on that.