I'm solving a system of two coupled partial integro-differential equations for two functions $ {\phi _0^\alpha (R,r')} $ and $ {\phi _0^\beta (R,r)} $: $$ \frac{{{\partial ^2}\phi _0^\alpha }}{{\partial {r^2}}} + f(R) \cdot \frac{{\partial \phi _0^\alpha }}{{\partial R}} + \frac{{{\partial ^2}\phi _0^\alpha }}{{\partial {R^2}}} + \left( {V\left( {r,R} \right) + \varepsilon _0^\alpha \left( R \right) + \int_a^b {v(r,r')|\phi _0^\beta (R,r'){|^2}dr'} } \right)\phi _0^\alpha = 0 $$ $$ \frac{{{\partial ^2}\phi _0^\beta }}{{\partial {r^2}}} + f(R) \cdot \frac{{\partial \phi _0^\beta }}{{\partial R}} + \frac{{{\partial ^2}\phi _0^\beta }}{{\partial {R^2}}} + \left( {V\left( {r,R} \right) + \varepsilon _0^\beta \left( R \right) + \int_a^b {v(r,r')|\phi _0^\alpha (R,r'){|^2}dr'} } \right)\phi _0^\beta = 0 $$ in which $f(R)$, $V(r,R)$ and $v(r,r')$ are known functions, $ {\varepsilon _0^\alpha \left( R \right)} $ and $ {\varepsilon _0^\beta \left( R \right)} $ are unknown, so it is also an eigenvalue problem.
The bounary conditions are Dirichlet boundary condition as: $$ \phi _0^\alpha \left( {A,r} \right) = {g_1}(r),\phi _0^\alpha (B,r) = {g_2}(r) $$ $$ \phi _0^\beta (R,a) = {h_1}(R),\phi _0^\beta (R,b) = {h_2}(R) $$
Since I'm a beginner in numerical solution of partial differential equations, I have no idea bout how to solve this PDE numerically. I have tried to use finite difference method but I don't know how to deal with the integral term and the eigenvalue term. Would Jacobi method or Gauss-Seidel method be helpful? Any hint is welcome. Thank you in advance for your help.