I wrote the code to do what you want. Please use it and test it. There are two slight differences:

1. Maple's range operator `..` is non-associative, and hence there's no way to write (a..b)..c without the parentheses. I don't think that it would be possible for Maplesoft to make it associative without affecting existing code because it works with lhs and rhs, just like `=`, and hence it inherently has two operands (either or both of which may be NULL).
2. My implementation requires inserting &.. between the rtable and the index specification if you want to use the new "skip" (i.e. "every nth") index specification that you requested. Unlike (1) above, this is not a "hard" requirement. It is possible to implement this new indexing without using any new operator between the rtable and the square brackets. But it is a little more complicated. One would need to ?overload the `?[]` operator.

Any standard-Maple index specification, no matter how complicated, should "pass through" my &.. operator and work the way it used to, so please test this aspect.

(*
 An implementation of "skip" or "every nth" rtable indexing: An attempt
 to implement some of the more-sophisticated Matlab-style indexing.
                                                                       *)
`&..`:= module()
option `Written 15-Apr-2013 by Carl Love.`;
local
     #Dimension bounds, set as each dimension is processed
    lo :: integer, hi :: integer,

     #Process the new type of index spec that this module introduces:
    #(a..b)..skip (also handles the case (..)..skip). Returns a list of the
    #actual indices to use.
    Skip:= proc(IndSpec :: range..posint)
        local a, b, k;
        a:= lhs(IndSpec);
        if a = (..) then (a,b):= (lo,hi) else (a,b):= op(a) end if;
        seq(k, k= `if`(a < 0 and lo = 1, hi+a+1, a)..`if`(b < 0 and lo = 1, hi+b+1, b), rhs(IndSpec))
    end proc,

    #Handle the dimension/index-spec pair for 1 dimension.
    DimProcess:= proc(Ind :: specfunc(anything, ``))
        local
            ind:= ``(op(2,Ind)),  #Index spec
            ret  #return value
        ;
        (lo,hi):= op(op(1,Ind));  #Bounds for this dimension
        ret:= expand(subsindets(ind, range..posint, Skip));
        if ind ::``(range..posint) then  [ret]  else  ret  end if     
    end proc,

     #Main code for &.. operator. The index spec for each dimension is paired
     #with that dimension's bounds.
    ModuleApply:= (A :: rtable, Inds :: list)->
        A[(DimProcess ~ (zip(``, [rtable_dims](A), Inds)))[]]
;
end module:     

IndexSkip.mw

The symbol e is not pre-defined in Maple for the purpose of input, although it does appear in Maple's output. You need to replace e^(-x) with exp(x) and replace e^(-10*x) with exp(-10*x). If you also solved the ODE with Maple, you'll need to resolve it using exp(-x) instead of e^(-x). Even though it may appear that Maple has given you a solution, that solution is just treating your input e as if it were any other variable.

Regarding your second plot, it is not clear what is meant by "plotting an equation", especially when that equation has only one variable. You could move the 300 to the left side of the equation, and then just plot the left side. The basic plot command does not handle equations.

Continuing right where your code left off:

ThePoints[1];
                           [0, 0, 0]
The first point, the origin, is not on the cylinder. I don't know why it is in the list. It messes up the matrix structure, and I am going to discard it for the rest of this work.

Group the points by coordinates in preparation for putting in matrices.
SortedPts:= sort(
   ThePoints[2..-1],
   (A,B)-> evalb(A[1] < B[1] or A[1] = B[1] and A[2] < B[2] or A[1..2] = B[1..2] and A[3] < B[3])
):
Convert from (x,y,z) form to (r, theta, z) form.
CylPts:= (P-> [sqrt(P[1]^2+P[2]^2), arctan(P[2],P[1]), P[3]]) ~ (SortedPts):
After the origin is discarded, the remaining points are naturally grouped 50 x 21 (21 values of z and 50 values of theta. r = 2 for all points, naturally). Check that r = 2 for all points:
remove(P->  fnormal(P[1]-2.0) = 0., CylPts);
                               []

Put into matrices
Z:= Matrix(21,50, (i,j)-> CylPts[(j-1)*21+i][3], datatype= float[8]):
Theta:= Matrix(21,50, (i,j)-> CylPts[(j-1)*21+i][2], datatype= float[8]):
R:= Matrix(21, 50, fill= 2., datatype= float[8]):

From the context of your question, I am guessing that you mean the base-10 logarithm. The symbol log in Maple, given without any other specification, means the base-e logarithm (also known as the natural logarithm and also known by the symbol ln). You can solve for x with this command:

solve(log10(x) = -4, x);

There is a command exactly for that. See ?fnormal.

You wrote:

And that works for small messages but what about a larger one like this: I see it stops after 6 "blocks".

sprintf(cat(WFO, cat(" ",WFO) $ words-1), sscanf(T, cat(WFI $ words))[]);

You are using the value of words that was set by my code for your previous "rainbow" string, which was 6. You need to re-execute this for your new string

# Count words
words:= iquo(length(S), wordlength, 'r'):
if r>0 then  words:= words+1  end if:

(using T instead of S), and you will get words=72. I used the variable words to count what you call "blocks" ("blocks" is a better name).

You wrote:

But do you mind explaining the WFI and WFO coding a little more?

I chose the names to stand for Word Format Input and Word Format Output. Could you ask a more specific question please? Also, please read ?cat, ?$ and the explanations of format codes at ?printf and ?scanf.

Using exact rational arithmetic, the size of the results is growing exponentially with each iteration, as the plot in the following worksheet shows. That means that even if the fastest algorithms that are theoretically possible are used, the time used for each iteration must grow at least exponentially. So, applying tricks like optimizing the code and multithreading might squeeze out another 1 or 2 iterations, but that's it.

What is your goal? Are you looking for an attractor or accumulation point? I think that you can find that with floating-point arthmetic. Can you do something with finite-field arithmetic? You could easily compute the 1000th iteration over 1000 different finite fields (with word-sized characteristics), but you'll never get the 1000th iteration in rational arithmetic (assuming the exponential trend continues).

Download ExpGrowth.mw

See the command ?CurveFitting,PolynomialInterpolation.

See the command ?LinearAlgebra:-GenerateMatrix. After that, you'll want to use ?LinearAlgebra,LinearSolve.

In addition to Adri's suggestions, you could set ?printlevel to a high value. For example

printlevel:= 100;

and then run your code.

Not having your code to work with, my first guess is that examplenr is not properly defined when you read it in another worksheet. I suggest that you construct the filename outside the plotsetup command and look at what it is. Something like this:

filename:= cat(filefolder, examplenr, ".gif");
plotsetup(gif, plotoutput= filename));

Anthony,

When looking at the description of the algorithm in, for example, Wikipedia, you have to think of y (the dependent variable) as a vector; and think of f (the function that evaluates the derivatives y') as a vector-valued function. Try the code below for your procedure Fourth. I kept the structure pretty close to what you had. If I were writing it from scratch, I would explicitly use Vectors.

Fourth:=proc(t,x,y,z,q,N)
local f,g,e,n,k1,k2,k3,k4,l2,l1,l3,l4,m1,m2,m3,m4,
   H, X, Y, Z
;
   f:= (x,y,z)-> a*x-q*x*y-nn*x*z;
   g:= (x,y,z)-> -b*y+k*x*y;
   e:= (x,y,z)-> -c*z+m*x*z;

   H:= h/2;

   for n from 0 to N do
      X:= x[n]; Y:= y[n]; Z:= z[n];
      k1:= f(X,Y,Z); l1:= g(X,Y,Z); m1:= e(X,Y,Z);

      X:= x[n]+H*k1; Y:= y[n]+H*l1; Z:= z[n]+H*m1;
      k2:= f(X,Y,Z); l2:= g(X,Y,Z); m2:= e(X,Y,Z);

      X:= x[n]+H*k2; Y:= y[n]+H*l2; Z:= z[n]+H*m2;
      k3:= f(X,Y,Z); l3:= g(X,Y,Z); m3:= e(X,Y,Z);

      X:= x[n]+h*k3; Y:= y[n]+h*l3; Z:= z[n]+h*m3;
      k4:= f(X,Y,Z); l4:= g(X,Y,Z); m4:= e(X,Y,Z);

      t[n+1]:= t[n] + h;
      x[n+1]:= x[n] + h/6*(k1 + 2*k2 + 2*k3 + k4);
      y[n+1]:= y[n] + h/6*(l1 + 2*l2 + 2*l3 + l4);
      z[n+1]:= z[n] + h/6*(m1 + 2*m2 + 2*m3 + m4)
   od

end proc:

We can use sscanf to break the input string into word-sized chunks, and then pass those words to sprintf so that they can be reassembled with the spaces between words. The syntax of sprintf and printf is identical. The reason that I chose sprintf is that there's no way to programmatically access the output of printf; it produces output for viewing only. 

restart;
S:= "IDIDNTSEETHATDOUBLERAINBOW":
wordlength:= 5:
# Construct scan format code
WFI:= sprintf("%%%ds", wordlength);

                         "%5s"

# Construct print format code
WFO:= sprintf("%%0.%ds", wordlength);

                          "%0.5s"

# Count words
words:= iquo(length(S), wordlength, 'r'):
if r>0 then  words:= words+1  end if:

sprintf(cat(WFO, cat(" ",WFO) $ words-1), sscanf(S, cat(WFI $ words))[]);
                            
               "IDIDN TSEET HATDO UBLER AINBO W"

Is this what you were looking for?

The error message about "points with fewer or more than 2 components" clearly provides a clue to the problem. The plot command is two-dimensional (2D), and thus cannot handle the 3D points that you try to give it in the plots axyz and ixyz. For these two plots, change the command plot to pointplot3d.

There's another bug in your code, more subtle and insidious. You use variable n several ways, some of which are in conflict. You use it as a parameter in the first differential equation de1, and as such it also appears in the f:= (x,y,z)-> ... in procedures ESys and Fourth. But you also use it as a local variable (for a loop index) in those procedures. You need to change n to something else, say nn, in four places: the definition of de1, the line where you set the numeric values of the parameters, and the two f:= (x,y,z)-> ... lines.

Consider adding some spaces to your code to make it easier to read.

If you have further questions about this piece of code, please upload a worksheet. This piece is a bit too long for an easy cut-and-paste by the reader (me in this case). After I cut-and-paste, I still have to manually block off the comments and remove the extra >s.

10^5 is not too many points for the plot renderer, but 10^6 is unwieldy. Execute these commands to create the plot:

nn:= nextprime(10^7):
zz:= 1:
N:= 10^5:
Z:= Matrix(N, 2):
for k to N do
   zz:= (zz^2+1) mod nn;
   Z[k,1]:= k:  Z[k,2]:= zz
end do:
plot(Z, style= point);

The plot will appear as a solid block of points. Put the mouse pointer in the plot and right click. Select Manipulator from the context menu, then Scale from the submenu. Select a rectangle on the plot which is horizontally narrow but vertically covers the full extent of the plot. Do this a few more times, until you're satisfied with the narrowness of the horizontal range. Then go to the manipuator submenu again and select Pan. Now the plot is horizontally scrollable.