I would do it this way.

> Phi:= (f(x*y) - g(x*y))*x*y;

> M:= f(x*y)*y; N:= g(x*y)*x;

> diff(M/Phi, y) - diff(N/Phi, x);

D(f)(x*y)/(f(x*y)-g(x*y))-f(x*y)/(f(x*y)-g(x*y))^2/x*(D(f)(x*y)*x-D(g)(x*y)*x)-D(g)(x*y)/(f(x*y)-g(x*y))+g(x*y)/(f(x*y)-g(x*y))^2/y*(D(f)(x*y)*y-D(g)(x*y)*y)

Hmm. The way Maple Primes displays this formula is not very nice.   Why does it switch D(f) to (f') (with a ' that is just one pixel on my screen)?  That should be

     D(f)(x y)      f(x y) (D(f)(x y) x - D(g)(x y) x)
  --------------- - ----------------------------------
  f(x y) - g(x y)                           2
                           (f(x y) - g(x y))  x

              D(g)(x y)      g(x y) (D(f)(x y) y - D(g)(x y) y)
         - --------------- + ----------------------------------
           f(x y) - g(x y)                           2
                                    (f(x y) - g(x y))  y

> normal(%);

  0

There isn't really a command for termwise integration.  You could do something like this.

> 1/sqrt(1-z^2) = convert(1/sqrt(1-z^2),FormalPowerSeries,z);

1/((1-z^2)^(1/2)) = Sum((2*k)!*4^(-k)/k!^2*z^(2*k),k = 0 .. infinity)

> int(lhs(%),z=0..x) = applyop(int,1,rhs(%),z=0..x);

arcsin(x) = Sum(int((2*k)!*4^(-k)/k!^2*z^(2*k),z = 0 .. x),k = 0 .. infinity)

> simplify(%);

arcsin(x) = Sum(x^(2*k+1)/(2*k+1)*(2*k)!*4^(-k)/k!^2,k = 0 .. infinity)

Are you looking for something like this?

> P:= map(t -> [t,1], select(isprime, [$2..100])):
  plots[pointplot](P, symbol=circle);

1)  See the section "THE SOLVING METHODS" in the help page ?dsolve,details.

2) Diff(M(x,y),y) - Diff(N(x,y),x) = 1, not 2

3) Are you talking about output or input?

For the numerical values of the roots in Maple, you could use fsolve.  Thus for all four roots at a = 2:

> fsolve(x^4 - 2*x - 1, x, complex);

    -.4746266176, -.4603551885-1.139317680*I, -.4603551885+1.139317680*I, 1.395336994

 Something else that's neat is rootlocus in the plots package.  Here are the paths in the complex plane of the four roots of x^4 - a*x - 1 as a goes from 0 to 10:

> plots[rootlocus](-x/(x^4-1),x, 0..10, scaling = constrained);

The expression you gave us has no variable z in it.  Did you mean to collect zg?

> collect(gg1, zg);

(2*tau[0]*D*a[1]*conjugate(q[2])^4*conjugate(ast)*b^2*A*c*conjugate(q[3])/(1+b*x[2*ast]+c*y[ast])^3-tau[0]*D*a[1]*conjugate(q[2])^2*conjugate(ast)*C*conjugate(q[3])^3*c^2/(1+b*x[2*ast]+c*y[ast])^3-tau[0]*D*a[1]*conjugate(q[2])^4*conjugate(ast)*C*conjugate(q[3])*b^2/(1+b*x[2*ast]+c*y[ast])^3-2*tau[0]*D*a[1]*conjugate(q[2])^3*conjugate(ast)*C*conjugate(q[3])^2*b*c/(1+b*x[2*ast]+c*y[ast])^3+tau[0]*D*a[1]*conjugate(q[2])^3*conjugate(ast)*b*A*c^2*conjugate(q[3])^2/(1+b*x[2*ast]+c*y[ast])^3+tau[0]*D*a[1]*conjugate(q[2])^5*conjugate(ast)*b^3*A/(1+b*x[2*ast]+c*y[ast])^3)*zg^4+(-tau[0]*D*a[1]*conjugate(q[2])^4*conjugate(ast)*b^2*A/(1+b*x[2*ast]+c*y[ast])^2-tau[0]*D*a[1]*conjugate(q[2])^2*conjugate(ast)*c*B*conjugate(q[3])^2*b/(1+b*x[2*ast]+c*y[ast])^2+tau[0]*D*a[1]*conjugate(q[2])^2*conjugate(ast)*C*conjugate(q[3])^2*c/(1+b*x[2*ast]+c*y[ast])^2-tau[0]*D*a[1]*conjugate(q[2])^3*conjugate(ast)*b*A*c*conjugate(q[3])/(1+b*x[2*ast]+c*y[ast])^2+tau[0]*D*a[1]*conjugate(q[2])^3*conjugate(ast)*C*conjugate(q[3])*b/(1+b*x[2*ast]+c*y[ast])^2)*zg^3+(-tau[0]*D*a*conjugate(q[2])^3*conjugate(ast)-tau[0]*D*a[1]*conjugate(q[2])^2*conjugate(ast)*C*conjugate(q[3])/(1+b*x[2*ast]+c*y[ast])+tau[0]*D*a[1]*conjugate(q[2])^3*conjugate(ast)*b*A/(1+b*x[2*ast]+c*y[ast]))*zg^2

If you want the result to be shown with powers of the variable zg increasing, you can use

> sort(%, zg, ascending);

Or for decreasing powers:

> sort(%, zg, descending);

I'm not sure what you mean by auto-rotation.  Do you mean something like this?

> with(plots): F:=plot3d(sin(x*y), x=-Pi..Pi, y=-Pi..Pi,colour=red):
 G:=plot3d(x + y, x=-Pi..Pi, y=-Pi..Pi,colour=green):
 H:=plot3d([2*sin(t)*cos(s), 2*cos(t)*cos(s), 2*sin(s)], s=0..Pi, t=-Pi..Pi):
 display({F, G, H},viewpoint="circleright");

Note: this requires Standard GUI.  The viewpoint option doesn't work in Classic.

Consider the polynomial P(x) = x^4 - a*x - 1 where a > 0.  Note that P'(x) = 4*x^3 - a < 0 for x < x0 = (a/4)^(1/3) and > 0 for x > x0, with P(x0) = -3/16*4^(2/3)*a^(4/3)-1 < 0.  Therefore there is exactly one real solution for x < x0 and exactly one for x > x0; the other two complex solutions must be a complex conjugate pair. 

Moreover, since P(0) < 0 < P(-1/a), one of the real solutions is in the interval (-1/a,0), and since P(1) < 0 < P(1 + a^(1/3)) the other real solution is in (1, 1 + a^(1/3)).

If you want to preserve spacing of Maple output that you paste in here, change the Format to "Formatted" instead of Normal.  For example, this is Normal:

> expand(P(a^(1/3)+1));
                               (2/3)      (1/3)
                      3 a + 6 a      + 4 a

while this is the same thing pasted in Formatted:

> expand(P(a^(1/3)+1));
                               (2/3)      (1/3)
                      3 a + 6 a      + 4 a

I might do it this way, taking advantage of the events option in dsolve(..., numeric) to stop when x or y hits 2 or x hits 0.

> xdot := diff(x(t),t) = x(t)-y(t):
  ydot := diff(y(t),t) = y(t)*(1-x(t)):
  interface(warnlevel=0):
  getarrows:= proc(x0, y0)
    local DEsol, tm,j;
    DEsol:= dsolve({xdot,ydot,x(0)=x0,y(0)=y0}, numeric,
      events=[[y(t)-2,halt],[x(t)-2,halt],[x(t),halt]],
      output=procedurelist);
    tm:= floor(subs(DEsol(9),t));
    seq(subs(DEsol(j), [x(t), y(t)]), j = 1 .. tm) 
    end proc:
  init:= map2(getarrows,0.2, [0.08, 0.085, 0.0873, 0.088]):
  DEtools[DEplot](
   {xdot,ydot}, 
   {x(t),y(t)}, 
   t=-20..0, 
   [[x(0)=0.9999,y(0)=0.9999],[x(0)=1.0001,y(0)=1.0001]], 
   scene=[x(t),y(t)], 
   x=0..2, y=0..2, 
   arrows=smalltwo, 
   dirfield=init, 
#   size=magnitude, 
#   color=magnitude,
   linecolor=black,
   thickness=6,
   numpoints=1000
 );

The straight line through point [x0, y0] with slope m has the equation y = y0 + m*(x - x0).

To get you started: what additional parameter do you need to know in order to find m_r(t)?

Once you have your equation, see the help page ?dsolve,numeric.

Yes, the coefficient of a^0 in the series is -1.  I don't know why convert(..., FormalPowerSeries,...) doesn't work with an indexed RootOf or other root selection mechanisms.  The positive root of x^4 - a*x - 1 (or at least the one that goes to +1 as a -> 0) should be

hypergeom([-1/12, 1/4, 7/12], [1/2, 3/4], -(27/256)*a^4)+(1/4)*a*hypergeom([1/6, 1/2, 5/6], [3/4, 5/4], -(27/256)*a^4)-(1/32)*a^2*hypergeom([5/12, 3/4, 13/12], [5/4, 3/2], -(27/256)*a^4)

if I'm not mistaken.

This all can be obtained from the Lagrange Inversion Theorem.  See e.g. groups.google.ca/group/sci.math.symbolic/msg/f278a2c4b127ce95

Maybe something like this?

> mkvec:= proc(s)
    uses StringTools;
    Vector(map(proc(c) if IsDigit(c) then parse(c) else convert(convert(c,name),`local`) end if end proc,
        Explode(s)))
  end proc;
  
 listplot([seq(i, i = 1 .. 50)], 
   tickmarks = [[seq(i = mkvec(StringTools[FormatTime]("%Y%b%d",timestamp=1260000000+60*60*24*i)), i = 1 .. 50,5)], 
      default], thickness = 3, color = green);
  

You can get the expression in radicals for the roots of the quartic if you set the environment variable

> _EnvExplicit:= true;

before calling solve.  However, it is not what I would call "reasonably simple".

> solve(v^4 - b*v - g = 0, v);

1/12*6^(1/2)*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)+1/12*(-(6*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-288*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*g-72*b*6^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3)/(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2))^(1/2),

1/12*6^(1/2)*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)-1/12*(-(6*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-288*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*g-72*b*6^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3)/(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2))^(1/2),

-1/12*6^(1/2)*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)+1/12*(-(6*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-288*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*g+72*b*6^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3)/(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2))^(1/2),

-1/12*6^(1/2)*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)-1/12*(-(6*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-288*(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2)*g+72*b*6^(1/2)*(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3)/(((108*b^2+12*(768*g^3+81*b^4)^(1/2))^(2/3)-48*g)/(108*b^2+12*(768*g^3+81*b^4)^(1/2))^(1/3))^(1/2))^(1/2)

 Maybe more interesting is the following.   First transform by scaling so the constant term is 1.

> Quartic:= v^4 - b*v - g;
   expand(eval(Quartic/g,v=u*g^(1/4)));

u^4-1/g^(3/4)*b*u-1

> convert(RootOf(%,u),FormalPowerSeries,b);

Sum(-pochhammer(7/12,k)*pochhammer(-1/12,k)*pochhammer(1/4,k)/pochhammer(3/4,k)/(2*k)!*(-27/64/g^3)^k*b^(4*k),k = 0 .. infinity)+Sum(1/4*1/g^(3/4)*pochhammer(5/6,k)*pochhammer(1/6,k)*(2*k)!^2*(-27/16/g^3)^k/k!^2/(4*k+1)!*b^(4*k+1),k = 0 .. infinity)+Sum(1/32*1/g^(3/2)*pochhammer(13/12,k)*pochhammer(5/12,k)*pochhammer(3/4,k)*(-27/64/g^3)^k/(2*k+1)!/pochhammer(5/4,k)*b^(4*k+2),k = 0 .. infinity)

 > g^(1/4)*value(%);

g^(1/4)*(-hypergeom([-1/12, 1/4, 7/12],[1/2, 3/4],-27/256*b^4/g^3)+1/4*1/g^(3/4)*b*hypergeom([1/6, 1/2, 5/6],[3/4, 5/4],-27/256*b^4/g^3)+1/32*1/g^(3/2)*b^2*hypergeom([5/12, 3/4, 13/12],[5/4, 3/2],-27/256*b^4/g^3))

To import the data, you might use ExcelTools[Import].  For example:

> M:= ExcelTools[Import]("c:/myFolder/myFile.xls", 1, "B5:C39");

To plot the data:

> with(plots): pointplot(M); P1:= %:

You must decide what type of function you want to use to describe the curve.  If you just want a straight line y = a*x + b, you could do this:
 

> X:= M[1..-1, 1]; Y:= M[1..-1, 2]:
  f:= Statistics[LinearFit]([1,x], X, Y, x);

Now to plot the best-fit line with the data:

> display([P1, plot(f, x = min(X) .. max(X))]);