This answer is, of course, assuming that you're talking about n -> infinity with k fixed.

This answer is, of course, assuming that you're talking about n -> infinity with k fixed.

Yes, Maple 9.5 (and only that release) has a bug that causes evalhf to be used even when Digits > evalhf(Digits).  A work-around is to make q2(x) an expression that can't be computed without using complex numbers (which evalhf in Maple 9.5 can't handle).

> q2:= (x) -> Re((x+I)*exp(x^2)*(1-erf(x))):
   Digits := 50:
   plot(q2(x), x = 0 .. 10);

Yes, Maple 9.5 (and only that release) has a bug that causes evalhf to be used even when Digits > evalhf(Digits).  A work-around is to make q2(x) an expression that can't be computed without using complex numbers (which evalhf in Maple 9.5 can't handle).

> q2:= (x) -> Re((x+I)*exp(x^2)*(1-erf(x))):
   Digits := 50:
   plot(q2(x), x = 0 .. 10);

Actually I should have made the initial conditions at t=0.

> ivp := { 0.01 * diff(y(t),t$2) - y(t)*diff(y(t),t) = 0, 
    y(0) = 0, D(y)(0) = -c } ;

> dsolve(ivp);

y(t) = 1/10*tan(5*t*(-c)^(1/2)*2^(1/2))*(-c)^(1/2)*2^(1/2)

> ysol:= simplify(convert(rhs(ysol),tanh)) assuming c > 0;

ysol := -1/10*tanh(5*t*c^(1/2)*2^(1/2))*c^(1/2)*2^(1/2)

> c0 := fsolve(eval(ysol,t=1) = -1, c = 0 .. infinity);

c0 := 50.00000000

> evalf(eval(ysol, {t=1, c = 50}), 50);

-.99999999999999999999999999999999999999999992559844

> evalf(eval(ysol, {t=0.9, c = 50}), 50);

-.99999999999999999999999999999999999999836119747516

> evalf(eval(ysol,{t=-1.1, c=50}), 50);

.99999999999999999999999999999999999999999999999658

Note that, the differential equation being autonomous, the solution values at t=0.9 and -1.1 are what we'd get at t=1 and -1 if we put the initial value at t=-0.1.

> plot([eval(ysol,c=50), eval(ysol,{t=t+0.1,c=50})],t=-1..1, 
     colour=[red,blue], axes=box);

The red curve is the solution with initial values y(0) = 0, y'(0) = 50.  The blue curve
is the solution with y(-0.1) = 0, y'(-0.1) = 50.  At t = 1 and t=-1 the y values differ from each other by less than 10^(-38).  Thus the boundary value problem specifying y values at -1 and 1 is extremely sensitive: very tiny changes in the boundary values produce big changes near t=0.  It's not at all surprising that numerical methods would have big troubles with this.

Actually I should have made the initial conditions at t=0.

> ivp := { 0.01 * diff(y(t),t$2) - y(t)*diff(y(t),t) = 0, 
    y(0) = 0, D(y)(0) = -c } ;

> dsolve(ivp);

y(t) = 1/10*tan(5*t*(-c)^(1/2)*2^(1/2))*(-c)^(1/2)*2^(1/2)

> ysol:= simplify(convert(rhs(ysol),tanh)) assuming c > 0;

ysol := -1/10*tanh(5*t*c^(1/2)*2^(1/2))*c^(1/2)*2^(1/2)

> c0 := fsolve(eval(ysol,t=1) = -1, c = 0 .. infinity);

c0 := 50.00000000

> evalf(eval(ysol, {t=1, c = 50}), 50);

-.99999999999999999999999999999999999999999992559844

> evalf(eval(ysol, {t=0.9, c = 50}), 50);

-.99999999999999999999999999999999999999836119747516

> evalf(eval(ysol,{t=-1.1, c=50}), 50);

.99999999999999999999999999999999999999999999999658

Note that, the differential equation being autonomous, the solution values at t=0.9 and -1.1 are what we'd get at t=1 and -1 if we put the initial value at t=-0.1.

> plot([eval(ysol,c=50), eval(ysol,{t=t+0.1,c=50})],t=-1..1, 
     colour=[red,blue], axes=box);

The red curve is the solution with initial values y(0) = 0, y'(0) = 50.  The blue curve
is the solution with y(-0.1) = 0, y'(-0.1) = 50.  At t = 1 and t=-1 the y values differ from each other by less than 10^(-38).  Thus the boundary value problem specifying y values at -1 and 1 is extremely sensitive: very tiny changes in the boundary values produce big changes near t=0.  It's not at all surprising that numerical methods would have big troubles with this.

Where did this x=0, x=1 in the first argument to plot come from?  It's not a syntax that I've ever seen before.  Is it documented somewhere?

The standard ways to produce the square would be

 plot([[0,0],[1,0],[1,1],[0,1],[0,0]], colour=black, axes=none);

or

 plots[display](plottools[rectangle]([0,0],[1,1], colour=black, style=line), axes=none);

By "shrink everything back until the hole is filled up again" do you mean transform (r, theta) to (r-a, theta) (in polar coordinates) for r >= a, where a is the radius of your hole?  If so, just plot max(r-a, 0) with coords=polar.

By "shrink everything back until the hole is filled up again" do you mean transform (r, theta) to (r-a, theta) (in polar coordinates) for r >= a, where a is the radius of your hole?  If so, just plot max(r-a, 0) with coords=polar.

Given a list R (which might be the result of Joe's code) consisting of small nonnegative integers, to plot a histogram you could use something like

Rmax := max (op(R));
counts:= [seq](numboccur(i,R), i = 0 .. Rmax);
Statistics:-ColumnGraph(counts, offset = - 0.4,    
  labels = ["Run length", "Frequency"]);

Given a list R (which might be the result of Joe's code) consisting of small nonnegative integers, to plot a histogram you could use something like

Rmax := max (op(R));
counts:= [seq](numboccur(i,R), i = 0 .. Rmax);
Statistics:-ColumnGraph(counts, offset = - 0.4,    
  labels = ["Run length", "Frequency"]);

Start with invlaplace in the inttrans package.  If you have a more specific question, we might be able to help.

Start with invlaplace in the inttrans package.  If you have a more specific question, we might be able to help.

Sometimes.  For example, consider this initial value problem:

> de := {(1+x)*(D@@2)(y)(x)+x*D(y)(x)+y(x)=0, y(0)=1, D(y)(0)=1};

Get the recurrence relation for the coefficients:

> gfun[diffeqtorec](de,y(x),a(n));

   {a(n)+n*a(1+n)+(n+2)*a(n+2),a(0)=1,a(1)=1}

Solve the recurrence:

> rsolve(%, a(n));

 (-1)^n*(1-3*n+n^2)/GAMMA(1+n)

So the series is

> y(x) = Sum(%*x^n, n=0..infinity);

Sometimes.  For example, consider this initial value problem:

> de := {(1+x)*(D@@2)(y)(x)+x*D(y)(x)+y(x)=0, y(0)=1, D(y)(0)=1};

Get the recurrence relation for the coefficients:

> gfun[diffeqtorec](de,y(x),a(n));

   {a(n)+n*a(1+n)+(n+2)*a(n+2),a(0)=1,a(1)=1}

Solve the recurrence:

> rsolve(%, a(n));

 (-1)^n*(1-3*n+n^2)/GAMMA(1+n)

So the series is

> y(x) = Sum(%*x^n, n=0..infinity);