rscalc[rm2m_] := rm2m + 2 m*(1 + Log[(rm2m)/(2 m)])
xmax = 700; xmin = -700; delx = 0.03; imax = Floor[(xmax - xmin)/delx];
rstar = xmax; m = 1.0; l = 2.0; rcour = 0.99; ntot = 3;
If[rstar > 4 m, r = rstar, r = 2 m*E^((rstar/(2 m)) - 1)];
rm2m = r - 2 m;
Do[rstest = rscalc[rm2m];
  drsdr = (rm2m + 2 m)/rm2m;
  rtempm2m = rm2m + (rstar - rstest)/drsdr;
  rm2m = rtempm2m;, {i, 7}];
mesh = NDSolve[{rm2mc'[rs] == rm2mc[rs]/(rm2mc[rs] + 2 m),
   rm2mc[rstar] == rm2m}, rm2mc, {rs, xmax, xmin}, AccuracyGoal -> 25,
   PrecisionGoal -> 25, WorkingPrecision -> 50]
Plot[Evaluate[rm2mc[lex] /. mesh], {lex, -20, 20}, PlotRange -> All]
xtest = -250;
b = Evaluate[rm2mc[xtest] /. mesh]
xchck = rscalc[b];
err = xchck - xtest
V[rm2m_] := (rm2m/(rm2m +
       2 m))*((l (l + 1))/(rm2m + 2 m)^2 + (2 m/(rm2m + 2 m)^3));
Plot[V[nax - 2 m], {nax, 0, 10}]
Vsx[x_] := V[Evaluate[rm2mc[x] /. mesh]];
Plot[Vsx[x], {x, -20, 20}]
XTab = Table[xmin + i*delx, {i, imax + 1}];