@want to be a permanent vegan I read your updated code, and I don't understand why odeplot is not satisfactory for you. Yes, you can use the output= Array(L) option to dsolve to pass in a list L of any number of independent-variable values. You'll receive as dsolve's return value an Array with the corresponding dependent-variables values (and some other stuff). Yes, you can make the spacing of the independent-variable values as fine as you want (for some extremely fine spacings you may need to adjust dsolve's error-control options). And yes, you can use these values to make a plot. But all these things can be handled better by odeplot, especially if you use its numpoints option.

@want to be a permanent vegan 

Please don't mix the MaplePrimes threads. That thread is about your if/piecewise problem. Please repost your latest Reply from that thread here, in this thread.

Neither the command-line interface (CLI) nor the print command is needed to realize this bug. Apparently, the bug is in the kernel. Any theories about how this arises? It seems that there's two distinct copies of t23.

@Muhammad Usman 

(vv: Thanks for explaining the raising to the 1/p power.)

Muhammad:

Please try this corrected translation of your algorithm. Note that it is divided into two procedures.

`mod/ReduceExponents`:= (f::polynom, p::posint)->
     #This procedure should only be used when all exponents of f are multiples of p.
     evalindets(f, `^`, t-> op(1,t)^iquo(op(2,t), p))
:

`mod/SFF`:= proc(f::polynom, p::posint)
local R:= 1, x:= indets(f, name)[], i, c, w, g, y, z;
     if x=NULL then return f mod p end if; #f is degree 0.
     g:= diff(f,x);
     if g=0 then return thisproc(ReduceExponents(f) mod p, p)^p end if;
     c:= Gcd(f,g) mod p;
     w:= Quo(f,c,x) mod p;
     for i while w <> 1 do
          y:= Gcd(w,c) mod p;
          z:= Quo(w,y,x) mod p;
          R:= R*z^i;
          w:= y;
          c:= Quo(c,y,x) mod p
     end do;
     R*thisproc(ReduceExponents(c) mod p, p)^p #returns R if c=1
end proc:

f:= x^11 + 2*x^9 + 2*x^8 + x^6 + x^5 + 2*x^3 + 2*x^2 + 1:
SFF(f) mod 3;

@Muhammad Usman 

Ah, yes, I didn't include the f^(1/p) and c^(1/p) steps. I don't know how to handle that. How can this be done without fully factoring the polynomial? I see that in the case you mention, you can just replace x^3 with x, but what about the general case?

@Just_Some_Guy 

In your worksheet, you have a vector of independent-variable values, X, which doesn't appear in your original Question. Also, what you called E---the independent variable in the original---has become dependent variable f in the above. Which is correct?

@Maple4evergr8 Before upgrading, you may want to confirm that that bug has indeed been fixed. I'm not sure. I'm sure that someone else reading this who has Maple 2015 can answer that question.

@Bendesarts You need to experiment with different initial conditions and different t-ranges until you find what you want. In the command below, the intial conditions are [1,1] and the t-range is 0..2.

DEtools[DEplot](EqSys, [x,z](t), t= 0..2, [[x,z](0)=~[1,1]]);

You need initial conditions: x(0) = ? y(0) = ?

@Preben Alsholm Thanks for the explanation. Although the work is done "by hand," the triangle explanation is sufficiently detailed to get a vote up from me. I know that Maple completely refuses to reverse the order of definite double integrals. I wonder if Mathematica or any other CAS has any capability at it for regions more complicated than rectangles.

@Preben Alsholm I have no doubt that your Answer is correct---the results can be verified with intsolve---but I wonder What's the rule of calculus that lets you change the limits of integration from x= 0..y to x= t..y? 

A variable p is used in the algorithm, in exponents only. I don't see where its value is set. So, what is p?

@J F Ogilvie Your reputation is next to your name on every post, next to your name at the top of the page (if you're logged in), and next to your name in many other places. Yours is currently 63.

@Preben Alsholm Thanks. Noted and corrected.

@Axel Vogt 

The links to a Mathematica marketing page and to a blog posted by Acer in a Reply to the Answer below suggest that Mathematica has implemented a general algorithm for highly oscillatory integrals rather than an ad hoc method to address the contest problem. I'm not a Mathematica user, so I can't test this.