plot(x^k, x=-1 .. 1, title = sprintf("The plot for value %d",k));

or even (when using Standard GUI)

plot(x^k, x=-1 .. 1, title = typeset("The plot of ",x^k));

It's interesting that you need that

 unprotect(D):

That is, the protection given to the global name D also affects the module export myD:-D.  I don't think that's a desirable feature.

It's interesting that you need that

 unprotect(D):

That is, the protection given to the global name D also affects the module export myD:-D.  I don't think that's a desirable feature.

for k in [ -1, -0.5, -0.25,0,25,0.5,1 ] do ...

for k in [ -1, -0.5, -0.25,0,25,0.5,1 ] do ...

In the Tools menu, choose Options, Security, check "Enable system/ssystem commands", and press "Apply Globally" (assuming you'll want changing directories to work all the time).  

But if you need the security of not allowing system commands (e.g. if you're going to be running code written by people you don't trust, who might do something nefarious with those commands), just include the path with the file name, e.g. in Doug's code

plotsetup( jpeg, plotoutput=`C:/mapleInputs/plot`||k||`.jpg`, etc.

In the Tools menu, choose Options, Security, check "Enable system/ssystem commands", and press "Apply Globally" (assuming you'll want changing directories to work all the time).  

But if you need the security of not allowing system commands (e.g. if you're going to be running code written by people you don't trust, who might do something nefarious with those commands), just include the path with the file name, e.g. in Doug's code

plotsetup( jpeg, plotoutput=`C:/mapleInputs/plot`||k||`.jpg`, etc.

It's a very bad idea to redefine Maple library functions, especially something like D that might be called by all sorts of other functions, which will now try to call your D instead. 
If you insist on doing it, you could try defining DD first, before changing D:

> DD:= eval(D):

But this is not going to work well.  In fact, it didn't work even without redefining D.

> DD(sin);

z -> DD(sin)(z)

which causes an infinite recursion if you try, e.g., DD(sin)(x).

Why do you want to make D into a number?  There's probably a better way...

It's a very bad idea to redefine Maple library functions, especially something like D that might be called by all sorts of other functions, which will now try to call your D instead. 
If you insist on doing it, you could try defining DD first, before changing D:

> DD:= eval(D):

But this is not going to work well.  In fact, it didn't work even without redefining D.

> DD(sin);

z -> DD(sin)(z)

which causes an infinite recursion if you try, e.g., DD(sin)(x).

Why do you want to make D into a number?  There's probably a better way...

I'm still not quite clear on what limit you are asking for.  You can't talk about a limit of cd/xd "for (x,c,q) constrained to the xd=0, cd=0, and qd=0 planes", because the denominator xd is 0 there.  You could ask for the limit of cd/xd to exist as, say, x -> x0 where f(x0) = c and f'(x0) = q.

Actually I think it's better to use t = sqrt(x) as the variable rather than x directly.  Otherwise there are problems with rationalizing expressions involving nested radicals.

> f:= x -> sqrt(x) - 1/10*x;
  solve(f(t^2)=c,t);

(25-10*c)^(1/2)+5, 5-(25-10*c)^(1/2), -5+(25-10*c)^(1/2), -(25-10*c)^(1/2)-5

We want a positive t.  This will require c < 2.5.  I'll choose the first solution.

> t0:= %[1];

> q0:= D(f)(t0^2) assuming c < 5/2;

q0 := 1/(2*((25-10*c)^(1/2)+5))-1/10

> indform := simplify((D(f)(t^2)-q0)/(f(t^2)-c)) assuming t>0;

indform := -5*((25-10*c)^(1/2)+5-t)/t/((25-10*c)^(1/2)+5)/(10*c+t^2-10*t)

> series(indform,t=t0,2);

5/2/((25-10*c)^(1/2)+5)^2/(25-10*c)^(1/2)+O((t-(25-10*c)^(1/2)-5))

 So the limit as x -> x0 = t0^2 with q = q0 fixed and c fixed is

> convert(%, polynom);

5/2/((25-10*c)^(1/2)+5)^2/(25-10*c)^(1/2)

The alleged solution is not a solution, because it makes the denominator x^0.5 - 0.1*x - c = 0.  In fact, you are explicitly requiring that denominator to be 0 because it is xd.  So how could there possibly be a solution?  Perhaps what you mean is this.

> sys:=map(numer,convert({xd,cd,qd},rational));

sys := {-(-5+x^(1/2)+10*q*x^(1/2))*c, 10*x^(1/2)-x-10*c, -50000*q^2*x+5000*x^(3/2)*q^2+2500*q*x-250*x^(3/2)*q-6250*x^(1/2)*q*c-30*x+3*x^(3/2)+155*x^(1/2)*c+25000*q*c-625*c}

> solve(%);

    {c = x^(1/2)-1/10*x, q = -1/10*(-5+x^(1/2))/x^(1/2), x = x}, {c = 0, q = q, x = 0}, {c = 0, q = q, x = 100}

Your problem with formatting is that your Input format has changed to Plain Text rather than Filtered HTML.  It sometimes does that for me too, for no apparent reason.  Click on Input format and you can change that.

For example, linear least squares could fit the curve y = a*x^2 + b*y + c where the unknown parameters are a,b,c, because this is a linear function of a,b and c, even though it is nonlinear in x.   But y = a*x^p + c would be nonlinear in the parameter p, so fitting it would be a nonlinear least-squares problem.

The result is 81 solutions, of which number 45 is equivalent to the single Maple 9.5 solution Alex quoted.  But number 63 is strange, because it has an extra variable that looks like O but seems to be an escaped local (it's not the global O):

{E = E, xi = RootOf(_Z^2+1), B[1] = B[1], B[2] = -(u[2]^2*B[1]-u[2]*B[1]*RootOf(_Z^2+1)*v[2]-2*u[1]*u[2]*B[1]+RootOf(_Z^2+1)*u[2]*E-RootOf(_Z^2+1)*u[1]*E)/u[1]/(u[2]+RootOf(_Z^2+1)*v[2]), p[1] = p[1], p[2] = -(-2*v[2]*u[2]^2*E^2*RootOf(_Z^2+1)+2*v[2]^3*RootOf(_Z^2+1)*u[2]^2*B[1]^2+4*u[1]^3*u[2]^2*O*B[1]^2-4*v[2]*u[2]^3*B[1]*E+2*u[2]^4*B[1]^2*O*u[1]-4*u[1]*B[1]^2*v[2]^2*u[2]^2-v[2]^3*RootOf(_Z^2+1)*p[1]*u[1]^2-6*u[2]^3*O*B[1]^2*u[1]^2-2*u[2]^3*O*B[1]^2*v[2]^2+2*v[2]*u[2]^4*B[1]^2*RootOf(_Z^2+1)+6*v[2]^3*RootOf(_Z^2+1)*u[1]^2*B[1]^2-2*v[2]^2*u[2]*E^2*O+2*v[2]^2*u[1]*E^2*O-4*v[2]^2*u[1]^3*O*B[1]^2-2*v[2]^2*u[2]*B[1]^2*u[1]^2+4*v[2]^3*E*u[1]*B[1]-v[2]^2*p[1]*u[1]^2*u[2]-6*v[2]^2*RootOf(_Z^2+1)*u[2]*O*E*B[1]*u[1]+6*v[2]^2*RootOf(_Z^2+1)*u[1]*u[2]*B[1]*E-6*v[2]^3*RootOf(_Z^2+1)*O*u[1]^2*B[1]^2-2*v[2]^3*RootOf(_Z^2+1)*O*u[2]^2*B[1]^2-6*v[2]^2*u[1]^2*B[1]*RootOf(_Z^2+1)*E+8*v[2]^3*RootOf(_Z^2+1)*O*u[2]*u[1]*B[1]^2-4*v[2]^3*B[1]*u[2]*E+6*v[2]^2*RootOf(_Z^2+1)*u[1]^2*E*O*B[1]+2*u[1]^2*u[2]^2*B[1]*RootOf(_Z^2+1)*E-2*v[2]^4*u[2]*O*B[1]^2-2*u[2]^4*u[1]*B[1]^2+6*u[2]^3*B[1]^2*u[1]^2-2*v[2]*RootOf(_Z^2+1)*u[1]*u[2]*E^2*O+8*v[2]*u[1]^3*u[2]*O*B[1]^2*RootOf(_Z^2+1)+2*v[2]^2*u[2]^3*B[1]^2-p[1]*u[1]^2*u[2]^3+2*v[2]*RootOf(_Z^2+1)*u[2]^2*E^2*O-8*v[2]*u[1]^3*u[2]*B[1]^2*RootOf(_Z^2+1)-8*v[2]*u[2]^3*B[1]^2*RootOf(_Z^2+1)*u[1]-4*v[2]^3*O*E*B[1]*u[1]+4*v[2]^3*u[2]*O*E*B[1]+2*v[2]^2*u[2]*O*B[1]^2*u[1]^2+4*u[2]^3*O*v[2]*E*B[1]-12*O*v[2]*E*B[1]*u[1]*u[2]^2+2*v[2]^4*u[1]*O*B[1]^2+2*v[2]*u[1]*E^2*RootOf(_Z^2+1)*u[2]+4*u[1]*O*B[1]^2*v[2]^2*u[2]^2+12*u[1]*B[1]*E*v[2]*u[2]^2-2*v[2]*u[2]^4*O*B[1]^2*RootOf(_Z^2+1)+14*v[2]*B[1]^2*u[1]^2*u[2]^2*RootOf(_Z^2+1)-v[2]*p[1]*u[1]^2*u[2]^2*RootOf(_Z^2+1)-4*u[1]^3*u[2]^2*B[1]^2-14*v[2]*u[1]^2*u[2]^2*O*B[1]^2*RootOf(_Z^2+1)+8*v[2]*u[2]^3*u[1]*O*B[1]^2*RootOf(_Z^2+1)+4*v[2]^2*u[1]^3*B[1]^2+2*v[2]^4*u[2]*B[1]^2-2*v[2]^4*B[1]^2*u[1]-2*v[2]^2*u[1]*E^2+2*v[2]^2*u[2]*E^2-2*RootOf(_Z^2+1)*u[1]*u[2]^3*B[1]*E-8*v[2]^3*RootOf(_Z^2+1)*u[2]*u[1]*B[1]^2-2*u[1]^2*u[2]^2*O*E*B[1]*RootOf(_Z^2+1)+2*u[2]^3*O*E*B[1]*u[1]*RootOf(_Z^2+1)+8*v[2]*u[1]^2*u[2]*O*E*B[1]-8*v[2]*u[1]^2*E*u[2]*B[1])/(RootOf(_Z^2+1)*v[2]^3+u[2]*v[2]^2+RootOf(_Z^2+1)*u[2]^2*v[2]+u[2]^3)/u[1]^2, rho[1] = -2*(3*u[1]*B[1]*E*v[2]-2*v[2]*u[2]*B[1]*E+u[2]^2*B[1]^2*O*u[1]+2*RootOf(_Z^2+1)*E*B[1]*v[2]^2+RootOf(_Z^2+1)*E^2*O*v[2]+RootOf(_Z^2+1)*O*B[1]^2*v[2]^3+2*RootOf(_Z^2+1)*u[1]^2*B[1]^2*v[2]+RootOf(_Z^2+1)*u[2]^2*B[1]^2*v[2]-2*u[2]*O*B[1]^2*u[1]^2-2*u[2]*O*B[1]^2*v[2]^2+3*u[1]*O*B[1]^2*v[2]^2-RootOf(_Z^2+1)*u[1]*u[2]*B[1]*E+2*u[2]*B[1]^2*u[1]^2-RootOf(_Z^2+1)*O*u[2]^2*B[1]^2*v[2]+2*RootOf(_Z^2+1)*O*u[2]*u[1]*B[1]^2*v[2]+RootOf(_Z^2+1)*u[2]*O*E*B[1]*u[1]-2*RootOf(_Z^2+1)*O*u[1]^2*B[1]^2*v[2]+2*u[2]*O*v[2]*E*B[1]-3*O*v[2]*E*B[1]*u[1]-2*RootOf(_Z^2+1)*u[2]*u[1]*B[1]^2*v[2]-2*RootOf(_Z^2+1)*E*O*B[1]*v[2]^2-RootOf(_Z^2+1)*E^2*v[2]-RootOf(_Z^2+1)*B[1]^2*v[2]^3-u[2]^2*u[1]*B[1]^2-3*u[1]*B[1]^2*v[2]^2+2*v[2]^2*u[2]*B[1]^2)/(-u[2]^3-u[2]*v[2]^2+RootOf(_Z^2+1)*v[2]^3+RootOf(_Z^2+1)*u[2]^2*v[2])/u[1]^2, rho[2] = -2*(3*u[1]*B[1]*E*v[2]-2*v[2]*u[2]*B[1]*E+u[2]^2*B[1]^2*O*u[1]+2*RootOf(_Z^2+1)*E*B[1]*v[2]^2+RootOf(_Z^2+1)*E^2*O*v[2]+RootOf(_Z^2+1)*O*B[1]^2*v[2]^3+2*RootOf(_Z^2+1)*u[1]^2*B[1]^2*v[2]+RootOf(_Z^2+1)*u[2]^2*B[1]^2*v[2]-2*u[2]*O*B[1]^2*u[1]^2-2*u[2]*O*B[1]^2*v[2]^2+3*u[1]*O*B[1]^2*v[2]^2-RootOf(_Z^2+1)*u[1]*u[2]*B[1]*E+2*u[2]*B[1]^2*u[1]^2-RootOf(_Z^2+1)*O*u[2]^2*B[1]^2*v[2]+2*RootOf(_Z^2+1)*O*u[2]*u[1]*B[1]^2*v[2]+RootOf(_Z^2+1)*u[2]*O*E*B[1]*u[1]-2*RootOf(_Z^2+1)*O*u[1]^2*B[1]^2*v[2]+2*u[2]*O*v[2]*E*B[1]-3*O*v[2]*E*B[1]*u[1]-2*RootOf(_Z^2+1)*u[2]*u[1]*B[1]^2*v[2]-2*RootOf(_Z^2+1)*E*O*B[1]*v[2]^2-RootOf(_Z^2+1)*E^2*v[2]-RootOf(_Z^2+1)*B[1]^2*v[2]^3-u[2]^2*u[1]*B[1]^2-3*u[1]*B[1]^2*v[2]^2+2*v[2]^2*u[2]*B[1]^2)/(-u[2]^3-u[2]*v[2]^2+RootOf(_Z^2+1)*v[2]^3+RootOf(_Z^2+1)*u[2]^2*v[2])/u[1]^2, u[1] = u[1], u[2] = u[2], v[1] = RootOf(_Z^2+1)*u[2]-u[1]*RootOf(_Z^2+1)+v[2], v[2] = v[2]}

I see no discontinuity near t = 0.038.  Can you upload your plot showing that?

I see no discontinuity near t = 0.038.  Can you upload your plot showing that?