You were there. Just expand about (x,y)=(1,-1) at truncation order 7 for total degree 6:

f:= (y^2) / (x^3) :
mtaylor(f,[x=1,y=-1],7);
                           2                              2
  2 - 3 x - 2 y + 6 (x - 1)  + 6 (y + 1) (x - 1) + (y + 1)

                    2                     3                     2
         - 3 (y + 1)  (x - 1) - 10 (x - 1)  - 12 (y + 1) (x - 1)

                     4                     3            2        2
         + 15 (x - 1)  + 20 (y + 1) (x - 1)  + 6 (y + 1)  (x - 1)

                     2        3             5                     4
         - 10 (y + 1)  (x - 1)  - 21 (x - 1)  - 30 (y + 1) (x - 1)

                     6             2        4                     5
         + 28 (x - 1)  + 15 (y + 1)  (x - 1)  + 42 (y + 1) (x - 1)

The exact sequence of statements that you are using is missing. In the following one, deduced from the fragments presented above, I do not see any such problem:

assume(  -1/2*R+M*G,  RealRange(-infinity,Open(0))):
((3*M*G*R-R^2)*R^2*R)*sqrt(-(3*M*G*R-R^2)^2*etc):
subs(R=1/U,%);
  (3*M*G/U-1/(U^2))/U^3*(-(3*M*G/U-1/(U^2))^2*etc)^(1/2)

If you do not write your system in Maple syntax, Maple won't give you any solution.

You can do subs(x=r,expression).

Where is this set of non-linear equations?

You may do the following:

J:=Int((exp(1/sqrt(2/delta)^2)/((t+delta*mu/(delta+1))^2*
(2*sqrt(2/delta)*sqrt(2*pi*delta*mu/(delta+1)))*t^(3/2))*
(t+delta*mu/(delta+1)))*exp(-(t*(delta+1)/(delta*mu)+delta*mu/
((delta+1)*t))/(2*sqrt(2/delta)^2)), t = 0 .. infinity);

subs(delta+1=a*delta*mu,J):
PDEtools:-dchange(t=u/a,%,[u]) assuming a>0:
simplify(%) assuming positive:
IntegrationTools:-Expand(%) assuming positive;

1/4*exp(delta)^(1/2)*delta^(1/2)*a^2/pi^(1/2)*Int(1/u^(3/2)*exp(-1/4*delta*u)*exp(-1/4*delta/u)/(u+1),u = 0 .. infinity)

So, you have to calculate the integral:

K:=Int(1/u^(3/2)*exp(-1/4*delta*u)*exp(-1/4*delta/u)/(u+1),
u = 0 .. infinity);

And it can be done by expanding 1/(1+u) in power series. Its n-th term is:

subs(1/(1+u)=(-1)^n*u^n,K):
Kn:=(simplify@value)(%) assuming n::posint;

Kn := 2*(-1)^n*BesselK(-1/2+n,1/2*delta)

and K becomes Sum(Kn,n=0..infinity).

use Statistics in Sample(RandomVariable(Normal(0, 1)),100)  end use;

in the GUI editor makes the job.

The rhs u^4-u+1 has no real root:

p := u^4-u+1:
[solve(p,u)]:
evalf~(%);

[.7271360845+.4300142883*I, -.7271360845+.9340992895*I, 
-.7271360845-.9340992895*I, .7271360845-.4300142883*I]

and is positive for real u. Hence, u(t) is monotonic increasing. And the time interval for u going from -infinity to +infinity is finite:

int(1/(u^4-u+1),u=-infinity..infinity):
Re(evalf(%));
                 2.683548244

So, depending on your initial condition, you may get a singularity "soon".

p:=sin(n*Pi/2):
simplify(subs(n=2*m,p)) assuming m::integer;
                                  0

simplify(subs(n=2*m+1,p)) assuming m::integer;

                                    m
                                (-1)

in the code X,Y,Min,Max.

If you write the ODE like diff(x(t),t)=f, you need to evaluate 'f'. So, e.g.:

f:=x(t)*sin(t):
dsys := {diff(x(t),t)=f,x(0)=1}: 
dsn := dsolve(dsys,numeric):
fn:=t->eval(f,dsn(t)):
fn(1);
                             1.332551265

The plot in Maxima 5.18 by:

plot2d(cos(2*x^2),[x,-10,10]);

(using gnuplot) is not so bad:

So, my suggestion for the calculation of the potential is to avoid the explicit expressions for Rp and Rm, which do not appear in your hand calculated result, but use instead function call forms, express  their derivatives back in terms of them, etc. Also, using the inert form of the functions may become convenient to keep under control the automatic simplification.

Indeed, the functionality of Maple is still rather primitive for this kind of manipulations, but something like these steps can be done:

vfX  := (X-1)*(2*Rp-Z)/(Rp^2*(Rp-Z)) - (X+1)*(2*Rm-Z)/(Rm^2*(Rm-Z));
vfX3:=subs(Rp=Rp(X,Y),Rm=Rm(X,Y),
X-1=diff(Rp(X,Y),X)*Rp(X,Y),X+1=diff(Rm(X,Y),X)*Rm(X,Y),vfX):
fp:=collect(diff(Rp(X,Y)*(Rp(X,Y)-Z),X),diff):
fm:=collect(diff(Rm(X,Y)*(Rm(X,Y)-Z),X),diff):
vfX4:=subs(op(1,fp)=%diff(Rp(X,Y)*(Rp(X,Y)-Z),X)/op(2,fp),
op(1,fm)=%diff(Rm(X,Y)*(Rm(X,Y)-Z),X)/op(2,fm),vfX3);
         d                              d
         -- (Rp(X, Y) (Rp(X, Y) - Z))   -- (Rm(X, Y) (Rm(X, Y) - Z))
         dX                             dX
 vfX4 := ---------------------------- - ----------------------------
           Rp(X, Y) (Rp(X, Y) - Z)        Rm(X, Y) (Rm(X, Y) - Z)

Method 'Contour' is not much developed. In particular:

int(cosh(a*x)/cosh(x),x=0..infinity,method=Contour);

returns unevaluated in Maple 13.

One way to compute this integral is:

J:=Int(cosh(a*x)/cosh(x),x=0..infinity):
convert(J,exp):
PDEtools:-dchange(x=ln(u),%): 
simplify(%) assuming 0<a,a<1;
value(%) assuming 0<a,a<1;

   1/2 LerchPhi(-1, 1, 1/2 - a/2) + 1/2 LerchPhi(-1, 1, a/2 + 1/2)