I'm solving a system of two coupled partial integro-differential equations as follow,
in which f(R), V(r,R) and v(r,r') are known functions, epsilon^alpha and epsilon^beta are unknown, so it is also an eigenvalue problem. The bounary conditions are known as phi^{alpha}(a,R), phi^{alpha}(b,R), phi^{alpha}(r,A) and phi^{alpha}(r,B) which can be express as four arrays, and phi^{beta} is the same(a, b, A, B are four constants of r and R respectively.)
I have tried to use method of lines and discretized the R dimension to transfer the PDEs to ODEs, but I found this method seems not work because MOL requires that the PDE problem is well-posed as an initial value problem in at least one dimension.
Since the PDE seems a little complicated, I have no idea about how to solve it numerically. Any hints will be helpful. Thank you in advance for your help.