I tried Ulrich Neumann's suggestion of the stackexchange link you gave and it works just fine.
To be explicit:
Vpot[p_, r_] :=
Piecewise[{{p[r]^2/(4*Pi),
p[r] < 217/500}, {-(87/250) + (4*p[r])/
5 + (87*(-1827 + 4400*p[r]))/(160000*Pi),
p[r] >= 109/
250}, {(86031603 + 640*Pi*(217 - 500*p[r])^2 -
405876800*p[r] + 500350000*p[r]^2 - 50000000*p[r]^3)/(800000*
Pi), True}}];
and
Upot[p_, r_] :=
Piecewise[{{1, p[r] <= 217/500}, {(957 + 320*Pi)/(200*p[r]),
p[r] >= 109/250}, {10007/4 + Pi*(800 - 1736/(5*p[r])) -
253673/(250*p[r]) - 375*p[r], True}}];
and using (for an example code snippet
...., Upot[p, r] --w^2/E^v[r])*p[r...,
within NDSolve's equation list. I got:
;
wprecision = 15
bosonstar[pcmax, 3] (*returns null in less than a couple seconds on a macbook laptop*)
but the value of rfin is defined
rfin (*returns 3.47146583566452*)