Numerically solving PDE system near singularity

Hello,
I’m trying to solve a PDE system near r=0, which is a singularity (please see attached file)
I chose to solve in the interval {r, \[Epsilon], R}, where \[Epsilon] is a positive cutoff.
The solution is very sensitive to changes in \[Epsilon], and for small enough values, the solution does not meet the boundary conditions.

Is there a better way to address this issue?

Thank you