I wonder whether there is a typo somewhere. I get a difference between the X component of the gradient of that hand calculated scalar potential:

sp1:=ln(Rp/Rm*(Rp-Z)/(Rm-Z));
                                  Rp (Rp - Z)
                        sp1 := ln(-----------)
                                  Rm (Rm - Z)

and vfX:

sp2:=subs(Rp=Rp(X,Y),Rm=Rm(X,Y),sp1):
s:=diff(Rp(X,Y),X)=(X-1)/Rp(X,Y),diff(Rm(X,Y),X)=
(X+1)/Rp(X,Y), diff(Rp(X,Y),Y)=Y/Rp(X,Y), 
diff(Rm(X,Y),Y)=Y/Rm(X,Y):
#Rp,Rm are introduced as function calls 
vf2X:=diff(sp2,X):
vf3X:=subs(s,vf2X):
fX:=op(indets(vf3X,specop(`*`,`+`))):
vf5X:=map(x->x*vf3X/fX,fX):
subs(X-1=Xp,X+1=Xm,vf5X):
collect(%,[Xp,Xm],simplify):
subs(Xp=X-1,Xm=X+1,Rp(X,Y)=Rp,Rm(X,Y)=Rm,%);

               (2 Rp - Z) (X - 1)   (2 Rm - Z) (X + 1)
               ------------------ - ------------------
                    2                 Rm Rp (Rm - Z)
                  Rp  (Rp - Z)

On the other hand, I recover vfY from the Y component of this gradient:

vf2Y:=diff(sp2,Y):
vf3Y:=subs(s,vf2Y):
fY:=op(indets(vf3Y,specop(`*`,`+`))):
vf5Y:=map(x->x*vf3Y/fY,fY):
vf6Y:=subs(Xp=X-1,Xm=X+1,Rp(X,Y)=Rp,Rm(X,Y)=Rm,%):
[selectremove(has,vf6Y,Rm)]:
collect(%,Y,simplify):
convert(%,`+`);

                      (2 Rm - Z) Y   (2 Rp - Z) Y
                    - ------------ + ------------
                        2              2
                      Rm  (Rm - Z)   Rp  (Rp - Z)

VF := (X,Y,u) -> VectorField( <Fx(X,Y,u), Fy(X,Y,u)>);
ScalarPotential( VF(X,Y,u) );
                                   2
                     - -------------------------
                         2            2    2 1/2
                       (X  - 2 X u + u  + Y )

It is a system of ODEs:

sol:=dsolve({eq1, eq2, BC1, BC2, BC3, BC4}, [u(y), b(y)]):

sol2:=convert(sol,trigh);
                     y                 sinh(30 y)
  sol2 := {b(y) = - ---- + _C4 + 1/900 ----------,
                     30                 cosh(30)

                             cosh(30 y)
        u(y) = 1/900 - 1/900 ----------}
                              cosh(30)

odetest(sol2,{eq1, eq2, BC1, BC2, BC3, BC4});

                                 {0}

Yes, there is a simpler way:

convert(sol3,sin,x):
collect(%,exp,simplify);

                   u(x, t) = 1/10 sin(x) exp(- t/2)

Presumably you have done:

K_op(psi);

                            2 m E psi
                            --------- = 0
                                  2
                              hbar

If so, at this point, 'psi' is an unassigned name, just as 'm' or 'E', representing something not depending on x and y. Hence, with zero derivatives. You can use as argument to 'K_op' a function call like psi(x,y), mainly for algebraic manipulations:

K_op(psi(x,y));
       / 2           \   / 2           \
       |d            |   |d            |   2 m E psi(x, y)
       |--- psi(x, y)| + |--- psi(x, y)| + --------------- = 0
       |  2          |   |  2          |            2
       \dx           /   \dy           /        hbar

or assign to 'psi' an specific expression for explicit computations:

psi:=sin(p*x)*sin(q*y):
K_op(psi):
solve(%,E);

                               2   2    2
                           hbar  (p  + q )
                           ---------------
                                 2 m

My guess.

As 'EllipticF' is the inverse of an elliptic function which is double-periodic, and x->0 corresponds to a pole of this elliptic function, which is "repeated" periodically in the complex plane, may be that 'MultiSeries:-limit' catches a pole shifted wrt that of the numerical methods.

If so, may be that some additional information should be given to 'MultiSeries:-limit' so that it catches the other pole.

You may do like this, starting with the third one:

data[3..-1];
                       [[2, 4], [3, 8], [4, 9]]

or also, in Maple 12:

data[3..];
                       [[2, 4], [3, 8], [4, 9]]

dq/dz=C*d/dz(dq/dz) is an ODE. Is this what you mean?

Proving deduced properties is the job of 'is'. But it is rather weak:

is(sin,continuous);
                                 true

is(cos,continuous);

                                 true

is(sin+cos,continuous);

                                 FAIL

Without a better prover, I would not expect too much.

Anyway, you may want to look at the subsection on Initial properties at the Maple Wiki.

issue occurs frequently. Eg:

Development tools

Emacs Mode

Vim mode

This package does not seem suitable for more than some simple computations in vector analysis. In particular 'DirectionalDiff' computes the directional derivative only for an scalar field.

Also, computations in any non-Cartesian system have severe problems.

but this thread seems relevant.

are you using? Could you upload your worksheet?

I would not say always, but frequently it proves convenient to "prepare" a complex expression before further computational steps. 

Prepare means, I think,  to transform to a mathematically equivalent expression of reduced computational complexity. May be using 'simplify', 'expand', it depends on the case. And it may require to use the assumptions.

You get close to g3 by:

(simplify@expand)(g1) assuming positive;


  /                                 /rho + 1\      / 1 \
  |                                 |-------|      |---|
  |                                 \  rho  /      \rho/
  \-inc[t + 1] R + a[t + 2] + a[t] R          delta

                   /rho + 1\      / 1 \\     /     / 1 \  / 1 \    \
                   |-------|      |---||     |     |---|  |---|    |
                   \  rho  /      \rho/|   / |     \rho/  \rho/    |
         + inc[t] R          delta     /  /  \delta      R      + R/
                                         /