Yes, evalc assumes all variables are real.  If you want s to be complex, substitute x+I*y for it.

Perhaps if you gave us a specific example of the type of system you're dealing with you might get some useful suggestions.

with(LinearAlgebra):
tmat1:=Matrix([[1,2,3],[2,4,5],[3,5,6]], shape=symmetric);
Eigenvectors(evalf(tmat1));

,

The third version, copied from the help page, used an "inert" version of differentiation.  Note that the d's are "greyed out".  The same effect can be obtained by typing Diff instead of diff.  To evaluate an expression containing inert operators, you can use the value command.

In the fourth version you apparently typed in d/dx.  Note that the d's are shown in italics.  This is not interpreted as differentiation, rather as a ratio of two variables named "d" and "dx".

I think you want to use the procedure option for dsolve.  See the help page ?dsolve,numeric,IVP.  Something like this:


> dsolve(numeric, number=1,
     procedure=proc(N,x,Y,YP) YP[1]:= evalf(Int(f(Y[1],p),p=p1..p2)) end proc,
     start= x0, initial = Array([y0]), procvars=[y(x)]);

Yes, there seems to be a bug here, perhaps in SolveTools:-UnwindRootOfs.  Assuming a > 0 doesn't help.  You might try the following work-around:

>  eq:=(a^(x+2))^(1/3)  =  (a^(x-5))^(1/2);
    expand(map(ln,eq)) assuming a>0,x::real;

> solve(%, x);

  19

The system of odes can be obtained as follows:

> R:= [r(t)*cos(theta(t)), r(t)*sin(theta(t))];
   A:= diff(R, t$2);
   des:= solve(expand(m*A + G*M*m*R/r(t)^3), 
             {diff(r(t),t$2), diff(theta(t),t$2)});

Now Kepler's second law says, essentially, that  , which is the rate at which the area is swept out, is constant.

> AM := r(t)^2*diff(theta(t),t);
   simplify(diff(AM,t), des);

           0

Converting to a list is quite inefficient if the Array is large.  You can try rtable_scanblock, e.g.

> a:= Array([11, 12, 13, 14, 15, 13, 16]);
   res:= table(); 
   rtable_scanblock(a, [],  
      proc(val,ind) global res; if val = 13 then  res[ind]:= NULL end if end proc);
   map(op,[indices(res)]);

Assuming that you have a renewal process, this is the Elementary Renewal Theorem.  See e.g. http://en.wikipedia.org/wiki/Renewal_theory or look up Renewal Theory in a probability text, e.g. Ross, "Introduction to Probability Models"

You only have residues at isolated singularities, not at branch points.

Let E = b^4 with b > 0.
In this case, I don't see how "any branch cut that is not on the real axis would be fine", but if you take your function as z* sqrt(b^2+z^2) * sqrt(1 - b^2/z^2) , with the principal branch of sqrt, you'll have branch cuts on the line segment [-b,b] of the real axis and [b,infinity)*I and (-infinity, -b]*I on the imaginary axis.  You could then take a contour such as the ellipse (x/a)^2 + (y/c)^2 = 1 with a > b > c and not meet a branch cut.

As far as I know, in all versions of Maple any expression you plot must evaluate to a numerical quantity after substituting a number for the variable.   You could "fake it" with something like this:

> plot(4*x, x=0..1, ytickmarks= [seq(i=a*i, i=0..4)], labels=[x,y]);

To obtain a square root function where the branch cut is at arg(z) = t, you can try

> mysqrt:= (t,z) -> sqrt(-z*exp(-I*t))/sqrt(-exp(-I*t));

Your system probably doesn't have closed-form solutions in general.  For particular values of the parameters and initial conditions, you can get numerical solutions using dsolve(..., numeric), and you can plot them using odeplot in the plots package.  For example:

> sys:= {diff(i(t),t) = alpha*B(t)*S(t)/(K+B(t))-r*i(t),
      diff(S(t),t) = -alpha*B(t)*S(t)/(K+B(t))+b*(N-S(t)),
      diff(B(t),t) = e*i(t)+g*B(t)};
   sys1:= eval(sys, {alpha=1, b=1, e=2,g=3, K = 4, N = 5, r=6});
   ic:= {i(0) = 1, S(0)=2, B(0)=0};
   sol:= dsolve(sys1 union ic, numeric);
   plots:-odeplot(sol,[[t,i(t)],[t,S(t)],[t,B(t)]],t=0..1, colour=[red,blue,green],
        legend=[i(t),S(t),B(t)]);

What exactly do you want to plot?  You can, for example, plot trajectories, something like this, for a system [x(k+1), y[k+1] =F([x(k), y(k)]):

>  F:= X -> [X[1] + (X[1]^2 - X[2]^2)/5, X[2] + X[1]*X[2]/5];
    Traj:= proc(X0, c, n)
       uses plots;
       local X, i;
       X[0]:= X0;
       for i from 1 to n do X[i]:= F(X[i-1]) end do:
       pointplot([seq(X[i],i=0..n)], colour=c);
     end proc;
    plots:-display(
     Traj~([[0.5,0],[0.5,0.5],[0,0.5],[-.5,.5],[-.5,0],[-.5,-.5],[0,-.5],[.5,-.5]],
          [red,orange,yellow,green,blue,cyan,magenta,black],10),
   axes=box, view=[-1..1, -1..1]);

  

It's sometimes very hard to predict how long something will take, especially since you haven't given us any details of what the equation and its solution look like.  But if it's taken over 15 hours without getting a result, I'd strongly suspect that either it's in an infinite loop or a combinatorial explosion that won't give you a result in our lifetimes.  Sometimes you can figure out what's going on by interrupting the calculation and using tracelast.   

You might get better results by using odetest.