1)  If small imaginary parts get smaller at higher values of Digits it is certainly an indication that there are roundoff issues.  It's not conclusive proof that the true value is real, but it would indicate that floating-point arithmetic is unlikely to decide that question.

2)  That's not the way I'd put it.  eval(J1, x=v) - eval(J1, x=0) is the integral from x=0 to v if J1 is an antiderivative of the integrand on the interval [0, v].   When that is not the case, you have to take into account the discontinuities in J1.   Maple's antiderivatives often do have jumps at some points where the integrand itself is smooth.  Definite integration by int(..., x = 0 .. v) will try to give a correct answer, taking this into account.  It doesn't always succeed.  On the other hand, evalf(Int(..., x = 0 .. v)) uses numerical integration methods, and usually is quite reliable.

>  with(Statistics):
   Coin1:= Sample(Bernoulli(4/5), 20000):
   Coin2:= Sample(Bernoulli(0), 20000):
   Coin3:= Sample(Bernoulli(1/2), 20000):
   A[0]:= 0: B[0]:= 0: C[0]:= 0:
   for t from 1 to 20000 do
     d1:= 2*round(Coin1[t])-1;
     d2:= 2*round(Coin2[t])-1;
     if A[t-1] mod 3 = 0 then A[t]:= A[t-1] + d2
     else A[t]:= A[t-1] + d1
     end if;
     if B[t-1] mod 3 = 1 then B[t]:= B[t-1] + d2
     else B[t]:= B[t-1] + d1
     end if;
     if C[t-1] mod 3 = Coin3[t] then
       C[t]:= C[t-1] + d2
     else
       C[t]:= C[t-1] + d1
     end if;
   end do:
  plot([[seq([i,A[i]],i=0..2000)],[seq([i,B[i]],i=0..2000)],[seq([i,C[i]],i=0..2000)]],
    symbol=point,style=point,colour=[red,blue,green]);

Note btw that in this case A[t] = B[t] + 2 for t >= 3.  Do you see why?

>  with(Statistics):
   Coin1:= Sample(Bernoulli(4/5), 20000):
   Coin2:= Sample(Bernoulli(0), 20000):
   Coin3:= Sample(Bernoulli(1/2), 20000):
   A[0]:= 0: B[0]:= 0: C[0]:= 0:
   for t from 1 to 20000 do
     d1:= 2*round(Coin1[t])-1;
     d2:= 2*round(Coin2[t])-1;
     if A[t-1] mod 3 = 0 then A[t]:= A[t-1] + d2
     else A[t]:= A[t-1] + d1
     end if;
     if B[t-1] mod 3 = 1 then B[t]:= B[t-1] + d2
     else B[t]:= B[t-1] + d1
     end if;
     if C[t-1] mod 3 = Coin3[t] then
       C[t]:= C[t-1] + d2
     else
       C[t]:= C[t-1] + d1
     end if;
   end do:
  plot([[seq([i,A[i]],i=0..2000)],[seq([i,B[i]],i=0..2000)],[seq([i,C[i]],i=0..2000)]],
    symbol=point,style=point,colour=[red,blue,green]);

Note btw that in this case A[t] = B[t] + 2 for t >= 3.  Do you see why?

To have an animation of a surface, you need three independent variables (two space and one time).  You have only two independent variables.  You could consider them both as space variables, producing a surface, or you can consider one as time and the other as space, producing an animation of a curve (as I did in my example).  

To have an animation of a surface, you need three independent variables (two space and one time).  You have only two independent variables.  You could consider them both as space variables, producing a surface, or you can consider one as time and the other as space, producing an animation of a curve (as I did in my example).  

The only difference is that you would use a*exp(b*x) instead of a*x^b.  Not a*e^(b*x), because Maple considers e to be an ordinary variable.

The only difference is that you would use a*exp(b*x) instead of a*x^b.  Not a*e^(b*x), because Maple considers e to be an ordinary variable.

If the lists have the same length, you can simply subtract:

> A-B;

 If not, you can use

> zip(`-`, A, B, 0);
 

which will pad the shorter one with 0's. 

If the lists have the same length, you can simply subtract:

> A-B;

 If not, you can use

> zip(`-`, A, B, 0);
 

which will pad the shorter one with 0's. 

indets finds the set of all sub-expressions of an expression that have a given type.  In this case the type is specfunc(anything,CURVES), which means any function call of the form CURVES(...).  Any plot command that produces curves will contain such objects.  Then op([1,1], ...) returns the first operand of the first operand of the set, i.e. the first operand of the CURVES function call, which is generally a list (but sometimes an Array) of points.

indets finds the set of all sub-expressions of an expression that have a given type.  In this case the type is specfunc(anything,CURVES), which means any function call of the form CURVES(...).  Any plot command that produces curves will contain such objects.  Then op([1,1], ...) returns the first operand of the first operand of the set, i.e. the first operand of the CURVES function call, which is generally a list (but sometimes an Array) of points.

You might be able to increase the stacklimit, if that is below the system's hard limit.  See ?kernelopts.
Alternatively, you might try using global variables rather than local variables of the procedures for storing your large data.

You might be able to increase the stacklimit, if that is below the system's hard limit.  See ?kernelopts.
Alternatively, you might try using global variables rather than local variables of the procedures for storing your large data.

The following procedure will generate an arrow of length L starting at the point on the solution Sol to the system Sys for t=T.  It is assumed that Sys is of the form {D(x)(t) = ..., D(y)(t) = ...} and Sol is the result of a call to dsolve({op(Sys), ICs}, numeric, output=procedurelist).

> arrowonsol:= proc(Sys, Sol, T, L)
         uses plots;
         local S, P, V;
         S:= Sol(T);
         P:= subs(S, [x(t), y(t)]);
         V:= evalf(subs(S, subs(Sys, [D(x)(t), D(y)(t)])));
         V:= L*V/sqrt(V[1]^2 + V[2]^2);
         plots[arrow](P,V,args[5..-1]);
    end proc;
  

For example:

>  with(plots):
   Sys:= [(D(x))(t) = -.1+x(t)^2-x(t)*y(t), (D(y))(t) = y(t)^2-x(t)^2-1];   
   Sol1:= dsolve([op(Sys), x(0)=0,y(0)=-1],numeric,output=procedurelist);
   Sol2:= dsolve([op(Sys), x(0) = 0.2, y(0) = -1.02], numeric, output=procedurelist);
   Curves:= odeplot(Sol1, [x(t),y(t)], t=-5..5),odeplot(Sol2, [x(t),y(t)], t=-5..5):
   display([Curves,
        seq(arrowonsol(Sys,Sol1,i,0.3,width=0.1),i=[-4,3]),
        seq(arrowonsol(Sys,Sol2,i,0.3,width=0.1),i=[-2.5,1.5])], 
     scaling=constrained);

The following procedure will generate an arrow of length L starting at the point on the solution Sol to the system Sys for t=T.  It is assumed that Sys is of the form {D(x)(t) = ..., D(y)(t) = ...} and Sol is the result of a call to dsolve({op(Sys), ICs}, numeric, output=procedurelist).

> arrowonsol:= proc(Sys, Sol, T, L)
         uses plots;
         local S, P, V;
         S:= Sol(T);
         P:= subs(S, [x(t), y(t)]);
         V:= evalf(subs(S, subs(Sys, [D(x)(t), D(y)(t)])));
         V:= L*V/sqrt(V[1]^2 + V[2]^2);
         plots[arrow](P,V,args[5..-1]);
    end proc;
  

For example:

>  with(plots):
   Sys:= [(D(x))(t) = -.1+x(t)^2-x(t)*y(t), (D(y))(t) = y(t)^2-x(t)^2-1];   
   Sol1:= dsolve([op(Sys), x(0)=0,y(0)=-1],numeric,output=procedurelist);
   Sol2:= dsolve([op(Sys), x(0) = 0.2, y(0) = -1.02], numeric, output=procedurelist);
   Curves:= odeplot(Sol1, [x(t),y(t)], t=-5..5),odeplot(Sol2, [x(t),y(t)], t=-5..5):
   display([Curves,
        seq(arrowonsol(Sys,Sol1,i,0.3,width=0.1),i=[-4,3]),
        seq(arrowonsol(Sys,Sol2,i,0.3,width=0.1),i=[-2.5,1.5])], 
     scaling=constrained);