@Markiyan Hirnyk I inverted the relationship between n and N from the original formulation. Thus I treat n rather than N as the symbolic variable, and I solve for it. This is handled by the lines

SolN:= [solve(b/c=n, N)]:
N:= SolN[1]:

I do this because it's easier to solve for an integer variable.

Your application depends on a code file, knot.txt, which hasn't been supplied.

@want to be a permanent vegan As we discussed earlier, that if-then-else clause needs to replaced by piecewise!

Can you find on the Web an example of the type of plot that you are talking about?

@Muhammad Usman You were very eager (to the point of being repetitive) to get an Answer, and now you've been silent for two days. So, does my most recently posted code work for you?

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