Using the solution to an equation in another function

I am trying to use the solution to $R^2 \left(1-\frac{1}{r}\text{Erf}[r/4]\right)=r^2$ for $r$ in the function

Potential$(R)=\frac{1}{R^2}\left(1-\frac{1}{r}\text{Erf}\left(r/4\right)\right)$

and plot this between $R=0$ and $R=5$. Essentially for each value of $R$, I want to solve the first equation for $r$ and then insert that $r$ into the second equation and plot. I use the code

T[R_] = First[r /. Solve[R^2 (1 - 1/r Erf[r/4]) == r^2, r, Reals]]
Potential[R_] = 1/R^2 (1 - 1/T[R] Erf[T[R]/4])
Plot[Potential[R], {R, 0, 5}]


However, T[R_] outputs $r$, which is not what I want, and then Potential[R_] outputs $(1 - (2 Erf[r/4])/r)/R^2$. Any ideas? I have tried using NSolve instead of Solve but this does not help.

• it returns r simply because Solve cant solve that. Do you imagine that equation should have an analytic solution? – george2079 Apr 15 '17 at 13:14
• I think it should, using NSolve then for example First[r /. NSolve[5^2 (1 - 1/r Erf[r/4]) == r^2, r, Reals]] gives -4.47771 – supercoolphysicist Apr 15 '17 at 13:51
• of course, NSolve can numerically handle many things, where Solve can not obtain an analytic soluton. – george2079 Apr 15 '17 at 16:40

Indeed, the plot can be obtained numerically.

T[R_?NumericQ] := First[r /. NSolve[R^2 (1 - 1/r Erf[r/4]) == r^2, r, Reals]]
Potential[R_?NumericQ] := 1/R^2 (1 - 1/T[R] Erf[T[R]/4])
Plot[Potential[R], {R, 0, 5}, PlotRange -> {0, 5}]


Note that Potential is singular at R == 0, so it might be more informative to plot.

Plot[R^2 Potential[R], {R, 0, 5}]


FYI you can do that plot parametrcially and avoid numerically solving for r :

RR[r_] = R /. Solve[R^2 (1 - 1/r Erf[r/4]) == r^2, R] // First
pot[r_] = 1/RR[r]^2 (1 - 1/r Erf[r/4])
ParametricPlot[{RR[r], RR[r]^2 pot[r]}, {r, -5, 0},
AspectRatio -> 1/GoldenRatio]


Or, another way could be find the solution for several numbers (x,y) and then interpolate as Interpolate will create a function for you.

interpF = Interpolation[datalist_as_x_y_pairs];


then

interpF[any_x] will get you the y value you are looking for