Both implicitplot and implicitplot3d, given an expression f, plot the curve or surface f = 0.  So (x^2 - y^2)/(x^2 + y^2) and x^2 - y^2 are the same from that point of view (except that (x^2 - y^2)/(x^2 + y^2) is not defined at [0,0]).  I repeat, what is "the limit" that you are trying to show does not exist?  Do you mean the limit of
(x^2-y^2)/(x^2+y^2)?  In that case, I think you want to look at something like

plot3d(  (x^2 - y^2)/(x^2 + y^2), x = -1 .. 1, y = -1 .. 1, axes = boxed);

What limit doesn't exist?  The intersection of the two surfaces z = x*y^2/(x^2+y^4) and x^2=y^2 consists of two curves: [x=t, y=t, z=t/(1+t^2)] and [x=t, y=-t, z=t/(1+t^2)] that behave quite nicely, intersecting at the origin (even though the limit of x*y^2/(x^2+y^4) as (x,y) -> (0,0) in two dimensions doesn't exist).

What limit doesn't exist?  The intersection of the two surfaces z = x*y^2/(x^2+y^4) and x^2=y^2 consists of two curves: [x=t, y=t, z=t/(1+t^2)] and [x=t, y=-t, z=t/(1+t^2)] that behave quite nicely, intersecting at the origin (even though the limit of x*y^2/(x^2+y^4) as (x,y) -> (0,0) in two dimensions doesn't exist).

Especially for something like implicitplot3d, this type of GUI zooming is not very useful, because the grid is not recalculated.  You just see bigger triangles. 

Especially for something like implicitplot3d, this type of GUI zooming is not very useful, because the grid is not recalculated.  You just see bigger triangles. 

> convert(1/(1+x^2),parfrac,x,I);

> A:= simplify(diff(%,x$n)) assuming n::posint;

A:=-1/2*I*GAMMA(1+n)*(-(-x+I)^(-1-n)+(-x-I)^(-1-n))

To get something like Gerald's form of it:

> evalc(%);

-1/2*GAMMA(1+n)*(exp(1/2*(-1-n)*ln(1+x^2))*sin((-1-n)*arctan(1,-x))-exp(1/2*(-1-n)*ln(1+x^2))*sin((-1-n)*arctan(-1,-x)))+1/2*I*GAMMA(1+n)*(exp(1/2*(-1-n)*ln(1+x^2))*cos((-1-n)*arctan(1,-x))-exp(1/2*(-1-n)*ln(1+x^2))*cos((-1-n)*arctan(-1,-x)))

But that can be simplified:

> simplify(%) assuming n::posint, x > 0;

GAMMA(1+n)*(1+x^2)^(-1/2-1/2*n)*sin((1+n)*(-arctan(1/x)+Pi))

However, better forms are available. 

> subs([I-x=1/rationalize(1/(I-x)),-I-x=1/rationalize(1/(-I-x))],A);

> simplify(convert(%,factorial)) assuming n::posint,  x > 0;

> simplify(%, size);

-1/2*I*n!*((x+I)*(-x-I)^n+(-x+I)*(-x+I)^n)*(1+x^2)^(-1-n)

> convert(1/(1+x^2),parfrac,x,I);

> A:= simplify(diff(%,x$n)) assuming n::posint;

A:=-1/2*I*GAMMA(1+n)*(-(-x+I)^(-1-n)+(-x-I)^(-1-n))

To get something like Gerald's form of it:

> evalc(%);

-1/2*GAMMA(1+n)*(exp(1/2*(-1-n)*ln(1+x^2))*sin((-1-n)*arctan(1,-x))-exp(1/2*(-1-n)*ln(1+x^2))*sin((-1-n)*arctan(-1,-x)))+1/2*I*GAMMA(1+n)*(exp(1/2*(-1-n)*ln(1+x^2))*cos((-1-n)*arctan(1,-x))-exp(1/2*(-1-n)*ln(1+x^2))*cos((-1-n)*arctan(-1,-x)))

But that can be simplified:

> simplify(%) assuming n::posint, x > 0;

GAMMA(1+n)*(1+x^2)^(-1/2-1/2*n)*sin((1+n)*(-arctan(1/x)+Pi))

However, better forms are available. 

> subs([I-x=1/rationalize(1/(I-x)),-I-x=1/rationalize(1/(-I-x))],A);

> simplify(convert(%,factorial)) assuming n::posint,  x > 0;

> simplify(%, size);

-1/2*I*n!*((x+I)*(-x-I)^n+(-x+I)*(-x+I)^n)*(1+x^2)^(-1-n)

Ummm, we don't know what w4_quiz_01 or w4_quiz_05 is.  You can get Fourier series coefficients using an integral.  Thus for a function f defined on the interval [0, 2*Pi] the Fourier coefficients (in the standard real form) are

a[n]=1/Pi*Int(f(x)*cos(n*x),x=0..2*Pi) for integers n >= 0

and

b[n]=1/Pi*int(f(x)*sin(n*x),x=0..2*Pi) for integers n >= 1

If doing this for symbolic n (rather than for particular integers) you'll probably want to make an assumption:

> assume(n::posint);

But watch out for exceptional values of n (in particular, the symbolic expression for a[n] often will not work for n=0).

It may be that fsolve will work better if you give it some help, in particular specify intervals for the variables.

> fsolve({.9999969940*y/((x+0.3006000000e-5)^2+y^2)^1.5
 +0.3006000000e-5*y/((x-.9999969940)^2+y^2)^1.5-y = 0.,
 .9999969940*(x+0.3006000000e-5)/((x+0.3006000000e-5)^2+y^2)^1.5
 +0.3006000000e-5*(x+.9999969940)/((x-.9999969940)^2+y^2)^1.5-x
 = 0.}, {x=-2..2,y=-1/2..1/2});

{x = 1.038151057, y = 0.}

It may be that fsolve will work better if you give it some help, in particular specify intervals for the variables.

> fsolve({.9999969940*y/((x+0.3006000000e-5)^2+y^2)^1.5
 +0.3006000000e-5*y/((x-.9999969940)^2+y^2)^1.5-y = 0.,
 .9999969940*(x+0.3006000000e-5)/((x+0.3006000000e-5)^2+y^2)^1.5
 +0.3006000000e-5*(x+.9999969940)/((x-.9999969940)^2+y^2)^1.5-x
 = 0.}, {x=-2..2,y=-1/2..1/2});

{x = 1.038151057, y = 0.}

That dsolve call produces an answer using an unevaluated integral.

t+Int(-1/(2*lambda*a*x1^2-2*lambda*a*x1+2*lambda*a*x2^2-2*lambda*a*x2-2*lambda*a*p*x1^2+2*lambda*a*p*x1-2*lambda*a*p*x2^2+2*lambda*a*p*x2-2*a*p*_a*x2-4*a*p*_a^2+4*a*p*_a-2*delta*_a),_a = `` .. x3(t))+_C1 = 0

For an explicit solution, try

> dsolve(eqn,x3(t),useint);

x3(t) = 1/4*(-p*a*x2+2*p*a-delta-tan(t*(-8*p*a^2*lambda*x1^2+8*p*a^2*lambda*x1-8*p*a^2*lambda*x2^2+8*p*a^2*lambda*x2+8*p^2*a^2*lambda*x1^2-8*p^2*a^2*lambda*x1+8*p^2*a^2*lambda*x2^2-8*p^2*a^2*lambda*x2-p^2*a^2*x2^2+4*p^2*a^2*x2-2*p*a*x2*delta-4*p^2*a^2+4*p*a*delta-delta^2)^(1/2)+_C1*(-8*p*a^2*lambda*x1^2+8*p*a^2*lambda*x1-8*p*a^2*lambda*x2^2+8*p*a^2*lambda*x2+8*p^2*a^2*lambda*x1^2-8*p^2*a^2*lambda*x1+8*p^2*a^2*lambda*x2^2-8*p^2*a^2*lambda*x2-p^2*a^2*x2^2+4*p^2*a^2*x2-2*p*a*x2*delta-4*p^2*a^2+4*p*a*delta-delta^2)^(1/2))*(-8*p*a^2*lambda*x1^2+8*p*a^2*lambda*x1-8*p*a^2*lambda*x2^2+8*p*a^2*lambda*x2+8*p^2*a^2*lambda*x1^2-8*p^2*a^2*lambda*x1+8*p^2*a^2*lambda*x2^2-8*p^2*a^2*lambda*x2-p^2*a^2*x2^2+4*p^2*a^2*x2-2*p*a*x2*delta-4*p^2*a^2+4*p*a*delta-delta^2)^(1/2))/p/a

> simplify(%, size);

x3(t) = 1/4*(-tan((8*p*((lambda*(x1^2-x1+x2^2-x2)-1/8*(x2-2)^2)*p-lambda*(x1^2-x1+x2^2-x2))*a^2-2*p*delta*(x2-2)*a-delta^2)^(1/2)*(t+_C1))*(8*p*((lambda*(x1^2-x1+x2^2-x2)-1/8*(x2-2)^2)*p-lambda*(x1^2-x1+x2^2-x2))*a^2-2*p*delta*(x2-2)*a-delta^2)^(1/2)-p*(x2-2)*a-delta)/p/a

To test whether this really works:

> odetest(%, eqn);

       0

That dsolve call produces an answer using an unevaluated integral.

t+Int(-1/(2*lambda*a*x1^2-2*lambda*a*x1+2*lambda*a*x2^2-2*lambda*a*x2-2*lambda*a*p*x1^2+2*lambda*a*p*x1-2*lambda*a*p*x2^2+2*lambda*a*p*x2-2*a*p*_a*x2-4*a*p*_a^2+4*a*p*_a-2*delta*_a),_a = `` .. x3(t))+_C1 = 0

For an explicit solution, try

> dsolve(eqn,x3(t),useint);

x3(t) = 1/4*(-p*a*x2+2*p*a-delta-tan(t*(-8*p*a^2*lambda*x1^2+8*p*a^2*lambda*x1-8*p*a^2*lambda*x2^2+8*p*a^2*lambda*x2+8*p^2*a^2*lambda*x1^2-8*p^2*a^2*lambda*x1+8*p^2*a^2*lambda*x2^2-8*p^2*a^2*lambda*x2-p^2*a^2*x2^2+4*p^2*a^2*x2-2*p*a*x2*delta-4*p^2*a^2+4*p*a*delta-delta^2)^(1/2)+_C1*(-8*p*a^2*lambda*x1^2+8*p*a^2*lambda*x1-8*p*a^2*lambda*x2^2+8*p*a^2*lambda*x2+8*p^2*a^2*lambda*x1^2-8*p^2*a^2*lambda*x1+8*p^2*a^2*lambda*x2^2-8*p^2*a^2*lambda*x2-p^2*a^2*x2^2+4*p^2*a^2*x2-2*p*a*x2*delta-4*p^2*a^2+4*p*a*delta-delta^2)^(1/2))*(-8*p*a^2*lambda*x1^2+8*p*a^2*lambda*x1-8*p*a^2*lambda*x2^2+8*p*a^2*lambda*x2+8*p^2*a^2*lambda*x1^2-8*p^2*a^2*lambda*x1+8*p^2*a^2*lambda*x2^2-8*p^2*a^2*lambda*x2-p^2*a^2*x2^2+4*p^2*a^2*x2-2*p*a*x2*delta-4*p^2*a^2+4*p*a*delta-delta^2)^(1/2))/p/a

> simplify(%, size);

x3(t) = 1/4*(-tan((8*p*((lambda*(x1^2-x1+x2^2-x2)-1/8*(x2-2)^2)*p-lambda*(x1^2-x1+x2^2-x2))*a^2-2*p*delta*(x2-2)*a-delta^2)^(1/2)*(t+_C1))*(8*p*((lambda*(x1^2-x1+x2^2-x2)-1/8*(x2-2)^2)*p-lambda*(x1^2-x1+x2^2-x2))*a^2-2*p*delta*(x2-2)*a-delta^2)^(1/2)-p*(x2-2)*a-delta)/p/a

To test whether this really works:

> odetest(%, eqn);

       0

Maple isn't using Matlab, but some Maple users are.  The Matlab package allows you to call Matlab from Maple, or to translate Matlab code to Maple.  One of the commands in this package is fft.  It uses Matlab to perform a discrete Fourier transform.  Maple also has its own Fast Fourier Transform command: FourierTransform in the DiscreteTransforms package.

From looking at several examples, it seems that when plotting a piecewise expression of the form

piecewise(x < 0, a(x), x = 0, b, x > 0, c(x))

with discont=true the point [0,b] is included if and only if a(0) <> b(0).  Whether this is the same as b does not matter.  For example:

> plot(piecewise(x < 0, 1, x = 0, 2, x > 0, 3), x = -1 .. 1, discont=true, axes = box);

> plot(piecewise(x < 0, 3, x = 0, 2, x > 0, 3-2*x), x = -1 .. 1, discont=true, axes = box);

Of course, to make sure the point is included you can always use pointplot.

The function f(x) = (-1)^x * x^(1/x) = exp(I*Pi*x + ln(x)/x) is not integrable on [1,infinity): note that f(x) ~ exp(I*Pi*x) as x -> infinity.  What does have a limit as N -> infinity (for integers N) is
int(f(x), x = 1 .. 2*N) .  Similarly, sum(f(n),n=1..2*N) has a limit as N -> infinity.