Last time, I asked this question

but after discussing with some friends, I have given up finding closed-form solutions. Now I have a simpler question.Let $g_i: i=1,2$ be $C^2 =C^{2}(-\infty,\infty)$ functions such that $0\leq g_1(x)\leq C_1$ and $0\leq g_2(x)\leq e^{x}C_2$ for some FIXED constants $C_1,C_2$. If $f_i : i=1,2$ are $C^2$ functions, let $\lambda$ be a fixed positive constants then define

$Af_i(x)=(\mu -\frac{\sigma^2}{2})\frac{\partial f_i(x)}{\partial x}+\frac{\sigma^2}{2}\frac{\partial^2 f_i(x)}{\partial x^2} +\lambda(g_i(x)-f_i(x)) $.

I am interested in finding NUMERICAL SOLUTIONS( i.e finding $f_1, f_2$ numerically) of the following systems of equations

$\min \{ \rho f_1(x) - Af_1(x), f_1(x)- f_2(x) +e^x D_1\} =0$

$\min \{ \rho f_2(x) - Af_2(x), f_2(x)- f_1(x) - e^x D_2\} =0$

where $\mu,\sigma,\rho, D_1>D_2$ are fixed positive constants

It would be highly appreciated if some one could give me some ideas on finding a scheme to solve this problem or at least a reference to this type of problem