@Christopher2222 Adding animatefield=true to the DEplot command makes an animation for each of the 90 frames in the outer animation. I don't think that is what you intended.

@Christopher2222 DEplot accepts the option maxfun= ... as does dsolve/numeric. Try maxfun=10^6 or whatever you think will work. In which setup are you using animatefield, i.e. which code?

@asa12 animate in Maple versions 12, 13, and 14 (but no later ones) is using subs in places where it ought to use eval.
## I now see in a reply to Kitonum that you actually use Maple 15. There my code works.
## Be aware that if you have axes = none, then what you see initially is just white space.
## Since the following may be useful to someone using earlier versions, I shall leave it as is.
##
One solution is to redefine `plots/animate` as follows.

restart;
## The first line redefines animate:
`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):
## Now test it:
sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]];
plots:-animate(plot,[sols[1..round(i)],thickness=3],i=1..nops(sols),paraminfo=false);

I just tried in Maple 12, and it works for me. I used this device when Maple versions earlier than 15 were the current ones.
As an alternative you can use this two step method:

restart;
sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]]; 
p:=proc(i) plot(sols[1..round(i)],_rest) end proc;
## Test the procedure p:
p(4,color=blue,thickness=3);
## Now animate p:
plots:-animate(p,[i,thickness=3,axes=none],i=1..nops(sols),paraminfo=false);

Using display/insequence as Kitonum does is perfectly fine. My point was just that animate can be used, if you want to.

Here is a procedure that parametrizes a polygon defined by a list of points given as [x,y].
(Since first posting this I changed it slighly so it can be used in any dimension).

par:=proc(pts::listlist,t::name) local d,n,i,j,r,s,p,px,py;
  d:=nops(pts[1]);
  n:=nops(pts);
  for i from 1 to n do
    j:=modp(i,n);
    r[i]:=pts[i]*(1-s)+pts[j+1]*s
  end do;
  p:=piecewise(op(op~([seq([t<i,eval(r[i],s=t-(i-1))],i=1..n)])),undefined);
  seq(evalindets(p,list,L->L[i]),i=1..d)
end proc;
## Now using CS as defined above we can do: 
px,py:=par(CS(12),t): 
plots:-animate(plot,[[px,py,t=0..T],thickness=7,color="SkyBlue"],
T=0..nops(CS(12)),paraminfo=false,axes=none,frames=100);

Here is an example of a space curve:
 

r:=rand(1..5.0);
L:=[seq([r(),r(),r()],i=1..100)]:
px,py,pz:=par(L,t):
plots:-spacecurve([px,py,pz],t=0..100,thickness=5,numpoints=5000);
plots:-animate(plots:-spacecurve,[[px,py,pz],t=0..T,thickness=5,numpoints=5000],T=0..100,frames=500);

Since it is Christmas here is a version with stars:
 

restart;
CS:=n->[seq((3+(-1)^k)*[cos(Pi/n*k),sin(Pi/n*k)],k=0..2*n)];
N:=8:
plots:-display(Array([seq(plots:-animate(plot,[CS(n)[1..round(i)],thickness=3,axes=none],i=1..nops(CS(n)),paraminfo=false),n=3..N)]),scaling=constrained);

It seems that I cannot show a gif file of the whole array, but I can give you one of the animations.
Here is CS(5):
 

plots:-animate(plot,[CS(5)[1..round(i)],thickness=3,axes=none],i=1..nops(CS(5)),paraminfo=false);

@rrbaldino An indication of the need for an assumption is given by the output from
Norm(`velθ`, 2);
    

Clearly the absolute values used are superfluous if the variables are real, which is an indication that Maple is open to the possibility that the variables are complex (as Maple generally is).

Could you upload a worksheet. We have no way of knowing what Vector[column](%id = 230612588) refers to.

@tsunamiBTP I believe it came with Maple 13. Cannot check though, I don't have that release anymore.

1. Don't you mean initial value problem instead of boundary value problem? Runge-Kutta methods are meant for initial value problems.
2. Once you have written the RK6 code for a single equation it shouldn't be much of a problem to extend it to systems.

Using the defining sum from the help pages we get:
 

restart;
S:=Sum(z^k*(product(pochhammer(n[i], k), i = 1 .. p))/(factorial(k)*(product(pochhammer(d[j], k), j = 1 .. q))), k = 0 .. infinity);
n:=[1, -1, 1/2];
d:= [-12,-3];
p:=nops(n);
q:=nops(d);
S;
value(S);

So the answer for z = 1 would be 71/72, which agrees with the numerical evaluation of the sum:

evalf(eval(S,z=1));
evalf(71/72);
# Interestingly, using the inert form with evalf gives the same answer (i.e. evalf(71/72)):

Hypergeom([1, -1, 1/2], [-12,-3], 1);
evalf(%);
Finally, you can try replacing z=1 with z=1.0:

hypergeom([1, -1, 1/2], [-12,-3], 1.0);

Thus I'm convinced you have found a bug. I shall report it as an SCR (Software Change Request).

@John Fredsted I couldn't agree more.

@denny The reason it doesn't make any difference in your case is that the ordering in sets is lexorder.
So try this:

restart;
{x(t),xdot(t)}; #order kept because x comes before xd in a dictionary.
restart;
{xdot(t),x(t)}; # same order in output as above.
restart;
{xm(t),xdot(t)}; #order reversed because xd comes before xm in a dictionary.
{xdot(t),xm(t)}; # order kept
### So use a list, when you want a particular order.

To be sure that the plot comes out with x along the horizontal axis and xdot along the vertical axis, use a list instead of a set.
Thus use [x(t), xdot(t)].
To see that it matters you could try your existing command with x(t) and xdot reversed just to see if you have control:
DEplot(sys,{xdot(t),x(t)}, ... and the rest ); #xdot along the vertical axis
Compare with
DEplot(sys,[xdot(t),x(t)], ... and the rest ); #x along the vertical axis

I got to think that you may wonder why the message `[Length of output exceeds limit of 1000000]`  came up, since Maple's own procedures are by default only printed as a meager skeleton.
To see why, try the simple code below.

restart;
ode:=diff(x(t),t)=x(t);
dsolve({ode,x(0)=1},numeric); # Only a skeleton is printed
interface(verboseproc); # default value 1
interface(verboseproc=2); # Setting the value to 2
dsolve({ode,x(0)=1},numeric); # Lots of stuff even for this simple problem.

The last value in HA, i.e. 0.8 can be handled by raising abserr and maxmesh some.
The following worked for me with the default setting of Digits at 10.

res[5] := dsolve(eval({Eq1, Eq2, Eq3, Eq4, IC1, IC2}, params union {K1 = lambda*HA[5]}), numeric,maxmesh=2048,abserr=1e-4,continuation=lambda);

It took a while and the message
`[Length of output exceeds limit of 1000000]`
came up. That is not an error message, so just ignore it.
The graphs of F and G at that value of K1:
plots:-odeplot(res[5],[[eta,F(eta)],[eta,G(eta)]]);