Another possible approach:
Given any value of E, you can solve your differential equation to find the corresponding \phi[E]. Conversely, given any \phi(r), you can perform the integral to compute E[\phi]. The combination of these is the fixed-point equation:
E[\phi[\mathcal{E}]] = \mathcal{E}
which is a nonlinear equation in a single variable that can be solved easily by a variety of root-finding algorithms.
The advantage of this is that you don’t need to implement an integro-differential solver directly — you just need an integration code for E[\phi] and a linear-ODE boundary-value solver for \phi[E], and then throw them at a root-finding algorithm.