I have no idea why your prof thinks using a "subroutine" to evaluate the polylogs is the key ingredient.
Subroutines are what Fortran had instead of procedures.  Maple  already has the capability of evaluating
polylogs.  A "Runge-Kutta like scheme" is used to solve an initial value problem.  Again, Maple can do
this with dsolve(..., numeric), with rather better methods than Runge-Kutta.  The real problem is that the
differential equation has no solutions satisfying the boundary condition y(infinity)=0.  All solutions will
go to infinity as r -> infinity.  I suggest telling that to your prof. 

I have no idea why your prof thinks using a "subroutine" to evaluate the polylogs is the key ingredient.
Subroutines are what Fortran had instead of procedures.  Maple  already has the capability of evaluating
polylogs.  A "Runge-Kutta like scheme" is used to solve an initial value problem.  Again, Maple can do
this with dsolve(..., numeric), with rather better methods than Runge-Kutta.  The real problem is that the
differential equation has no solutions satisfying the boundary condition y(infinity)=0.  All solutions will
go to infinity as r -> infinity.  I suggest telling that to your prof. 

Well, then you do want a 3d plot.  But the world should not be flat, it should be a globe.
I think it may be best to plot the continent outlines as curves rather than using an image.

Well, then you do want a 3d plot.  But the world should not be flat, it should be a globe.
I think it may be best to plot the continent outlines as curves rather than using an image.

I think what you want is

> mtaylor(f(x+dx,y+dy,t+dt),[dx,dy,dt],3);

f(x,y,t)+D[1](f)(x,y,t)*dx+D[2](f)(x,y,t)*dy+D[3](f)(x,y,t)*dt+1/2*D[1,1](f)(x,y,t)*dx^2+dx*D[1,2](f)(x,y,t)*dy+dx*D[1,3](f)(x,y,t)*dt+1/2*D[2,2](f)(x,y,t)*dy^2+dy*D[2,3](f)(x,y,t)*dt+1/2*D[3,3](f)(x,y,t)*dt^2

I think what you want is

> mtaylor(f(x+dx,y+dy,t+dt),[dx,dy,dt],3);

f(x,y,t)+D[1](f)(x,y,t)*dx+D[2](f)(x,y,t)*dy+D[3](f)(x,y,t)*dt+1/2*D[1,1](f)(x,y,t)*dx^2+dx*D[1,2](f)(x,y,t)*dy+dx*D[1,3](f)(x,y,t)*dt+1/2*D[2,2](f)(x,y,t)*dy^2+dy*D[2,3](f)(x,y,t)*dt+1/2*D[3,3](f)(x,y,t)*dt^2

And here's one way to get a colour image as a 2D plot.  Warning again: this produces a huge plot object.  That's probably why ImageTools:-Preview makes a 3d plot rather than 2d.

> FileLocation:= cat(kernelopts(datadir), "/images/tree.jpg");
  A:= ImageTools:-Rotate(ImageTools:-Read(FileLocation),90):
  for i from 1 to 3 do
    P[i]:= plots[listdensityplot](A[1..-1,1..-1,i],style=patchnogrid);
    C[i]:= op(indets(P[i],specfunc(anything,COLOR)));
    L[i]:=convert(convert(C[i],list)[2..-1],Vector) 
  end do:
  CC:= COLOR(RGB,seq(L[(i-1) mod 3+1][i], i = 1 .. nops(C[1])-1)):
  subs(C[1]=CC,P[1]);

And here's one way to get a colour image as a 2D plot.  Warning again: this produces a huge plot object.  That's probably why ImageTools:-Preview makes a 3d plot rather than 2d.

> FileLocation:= cat(kernelopts(datadir), "/images/tree.jpg");
  A:= ImageTools:-Rotate(ImageTools:-Read(FileLocation),90):
  for i from 1 to 3 do
    P[i]:= plots[listdensityplot](A[1..-1,1..-1,i],style=patchnogrid);
    C[i]:= op(indets(P[i],specfunc(anything,COLOR)));
    L[i]:=convert(convert(C[i],list)[2..-1],Vector) 
  end do:
  CC:= COLOR(RGB,seq(L[(i-1) mod 3+1][i], i = 1 .. nops(C[1])-1)):
  subs(C[1]=CC,P[1]);

Depending on what you mean by that, a "completely even set of points" on the sphere may not be possible.  For N points with a statistically uniform distribution on the sphere in R^n, try something like this.

> with(Statistics):
  Z:= RandomVariable(Normal(0,1));
  L:= Sample(Z, N*n);
  R:= [seq(add(L[i]^2,i=j*n+1 .. (j+1)*n)^(1/2),j = 0 .. N-1)]; 
  A:= <seq(L[j*n+1 .. (j+1)*n]/R[j+1], j=0..N-1)>;

The points are the rows of the resulting Matrix A.