As you write your problem in integral form then I understand the integration is not easily expressible as compact formula after it is performed. At least when the functions b(.) and \sigma(.,.) are not specified explicitly. I also assume that both unknowns, \varepsilon_m and \chi are just real numbers, not functions. If this is the case, then I would recommend to use interval methods to solve this problem numerically. (More on interval computing at https://en.wikipedia.org/wiki/Interval_arithmetic).
Without going too deeply into details: interval calculations make possible to evaluate a range of (real) expression when bounds of all its components are known. When you have an idea in what ranges your unknowns are expected to be found, the the problem boils down to systematic narrowing of those ranges until they are sufficiently narrow. The many existing software tools to handle interval problems are also mentioned at given Wikipedia address.
I hope you'll become a lover of interval methods once you try them.
I received an excellent algorithm from Prof Polyanin, but I thought of a theoritical explicit solution or a specific numerical method for the integral equations systems.
That's right, b(.) and sigma(.) are arbitrary, and the ranges of epsilon[0-1%] and chi[0.05-1].
I will read what it is written in Wikipedia, and if I have more questions, I won't hesitate