My code works in Maple 13 Standard, but not Classic: in Classic, pointplot  doesn't allow a list of colours.  This is a bug.  Leave out the colour option (or change to colour=red) and it will work.

In your example, you're getting one point per frame because that's what you gave pointplot: c[i] is just one point.   You might try

display(seq(pointplot(c[1..i]),i=1..nops(c)), insequence=true);

You might find this interesting: an animation of points running around the attractor.

> sys:= eval({ross_x,ross_y,ross_z, x(0)=1, y(0)=1,z(0)=1},params);
  with(plots):
  bg:= odeplot(dsolve(sys, numeric, maxfun=infinity),
    [x(t),y(t),z(t)],t=1000..2000, numpoints=15000):
  M:= dsolve(sys, numeric, maxfun=infinity, 
     output=Array([seq(i,i=1000 .. 1040, 0.1)]))[2,1];
  for i from 0 to 39 do
    frame[i]:= display([bg, 
       pointplot3d([seq(M[1+i+40*j,2..4],j=0..9)], 
          colour=black, symbol=solidsphere)])
  end do:
  display([seq(frame[i],i=0..39)],insequence=true);

It's not so clear to me that it's possible to "overcome" the singularity.  It appears that at approximately t=2.6991460, the derivatives of theta[1] and theta[2] become infinite.
For example,

> ans1(2.699145);

[t = 2.699145, theta[1](t) = 158.260352315183696, diff(theta[1](t), t) = -945782.489321760136, theta[2](t) = 122.244130903344143, diff(theta[2](t), t) = -945782.152519616532]

You could try using a "classical" method with fixed step size: then your solution will not go to infinity at a finite time, but it will still become extremely large.  For example:

> RK:=dsolve({eq1, eq2, theta[1](0) = 170.8, theta[2](0) = 124.1,
 (D(theta[1]))(0) = 0, (D(theta[2]))(0) = 0}, 
 [theta[1](t), theta[2](t)], numeric, 
 method=classical[rk4], stepsize=0.001);

> RK(3);

[t = 3., theta[1](t) = -7.91552425031953828*10^65, diff(theta[1](t), t) = -2.64291293833512376*10^66, theta[2](t) = -7.91552425031953828*10^65, diff(theta[2](t), t) = -2.64291293833512376*10^66]

In case you missed it, those are numbers on the order of 10^67.  I very much doubt that these are physically meaningful.

I think you should re-examine your model to see if the differential equation is correct, and if so what is happening to make it "blow up", and whether there are other variables that could be used that would not blow up.

It's not so clear to me that it's possible to "overcome" the singularity.  It appears that at approximately t=2.6991460, the derivatives of theta[1] and theta[2] become infinite.
For example,

> ans1(2.699145);

[t = 2.699145, theta[1](t) = 158.260352315183696, diff(theta[1](t), t) = -945782.489321760136, theta[2](t) = 122.244130903344143, diff(theta[2](t), t) = -945782.152519616532]

You could try using a "classical" method with fixed step size: then your solution will not go to infinity at a finite time, but it will still become extremely large.  For example:

> RK:=dsolve({eq1, eq2, theta[1](0) = 170.8, theta[2](0) = 124.1,
 (D(theta[1]))(0) = 0, (D(theta[2]))(0) = 0}, 
 [theta[1](t), theta[2](t)], numeric, 
 method=classical[rk4], stepsize=0.001);

> RK(3);

[t = 3., theta[1](t) = -7.91552425031953828*10^65, diff(theta[1](t), t) = -2.64291293833512376*10^66, theta[2](t) = -7.91552425031953828*10^65, diff(theta[2](t), t) = -2.64291293833512376*10^66]

In case you missed it, those are numbers on the order of 10^67.  I very much doubt that these are physically meaningful.

I think you should re-examine your model to see if the differential equation is correct, and if so what is happening to make it "blow up", and whether there are other variables that could be used that would not blow up.

I'm not sure what kind of animation you want.  Do you want a different random sample for each frame of the animation, or one random sample which is then followed for different frames?

Well, for example here's an animation of 100 particles undergoing Brownian motion.
For this kind of animation, it may be best to construct frames one by one, then put them together using display(..., insequence=true).  Thus:

> with(plots):
  with(Statistics):
  X:= Matrix(100,2,datatype=float[8]);
  colours:= [seq(COLOUR(HUE,0.01*i),i=1..100)]:
  frame[0]:= pointplot(X, colour=colours):
  for i from 1 to 60 do
    X:= X + 
       ArrayTools[Reshape](Sample(Normal(0,0.1),200),100,2);
    frame[i]:= pointplot(X,colour=colours)
  end do:
  display([seq(frame[i],i=0..60)], insequence=true);
 

I'm not sure what kind of animation you want.  Do you want a different random sample for each frame of the animation, or one random sample which is then followed for different frames?

Well, for example here's an animation of 100 particles undergoing Brownian motion.
For this kind of animation, it may be best to construct frames one by one, then put them together using display(..., insequence=true).  Thus:

> with(plots):
  with(Statistics):
  X:= Matrix(100,2,datatype=float[8]);
  colours:= [seq(COLOUR(HUE,0.01*i),i=1..100)]:
  frame[0]:= pointplot(X, colour=colours):
  for i from 1 to 60 do
    X:= X + 
       ArrayTools[Reshape](Sample(Normal(0,0.1),200),100,2);
    frame[i]:= pointplot(X,colour=colours)
  end do:
  display([seq(frame[i],i=0..60)], insequence=true);
 

What do you mean by "discontinuity checking" in implicitplot, if not what signchange does?

What do you mean by "discontinuity checking" in implicitplot, if not what signchange does?

There is no such limitation.  The only limitation on lists is that you can't assign to a member of a long list.  But there's no problem with constructing a long list using seq or map or zip. 

You could try e.g.

> f:= [seq(sin(i),i=1..1000)]:
  id:= t -> t;
  pointplot3d(map([id, id^2, rand(-50..50)], f));

There is no such limitation.  The only limitation on lists is that you can't assign to a member of a long list.  But there's no problem with constructing a long list using seq or map or zip. 

You could try e.g.

> f:= [seq(sin(i),i=1..1000)]:
  id:= t -> t;
  pointplot3d(map([id, id^2, rand(-50..50)], f));

What implicitplot does is detect changes in sign.  For the equation y = tan(x), the difference between the sides changes sign as x crosses an odd multiple of Pi/2.  That's what produces the vertical lines.  If you only want sign changes where the value goes through 0, you can use the option signchange=false.

> plots[implicitplot](y=tan(x),x=-5..5,y=-5..5,
    signchange=false,gridrefine=2);

What implicitplot does is detect changes in sign.  For the equation y = tan(x), the difference between the sides changes sign as x crosses an odd multiple of Pi/2.  That's what produces the vertical lines.  If you only want sign changes where the value goes through 0, you can use the option signchange=false.

> plots[implicitplot](y=tan(x),x=-5..5,y=-5..5,
    signchange=false,gridrefine=2);

If you don't predefine A and B, what you get are tables, which can't be added directly.  You could try

> C:= Array(0..2,0..2,A) + Array(0..2,0..2,B);

If you don't predefine A and B, what you get are tables, which can't be added directly.  You could try

> C:= Array(0..2,0..2,A) + Array(0..2,0..2,B);

I presume your web browser has been set up to block popup windows.  You should place www.mapleprimes.com on the list of sites from which popup windows are allowed.  If you're using Firefox, choose Tools, Options..., Content, click on the Exceptions... button to the right of "
"Block popup windows", and www.mapleprimes.com and click Allow.