Solving this system of ODE

I have the following system of equations

$$\rho V(u, \epsilon^i) = F(u, \epsilon^i) + V_u(u, \epsilon^i)g(u, V(u, \epsilon^i) + \lambda^i \left(V(u, \epsilon^{-i}) - V(u, \epsilon^i)\right)$$

with $i \in \{0,1\}$, so that $\epsilon$ can take two values and $g(u, V)$ known and initial $u(0), \epsilon$ given. $g(u,V)$ is the law of motion for $u$, $\dot u(u, \epsilon^i)$ if you will - but it depends on the value of $V$. This looks like a HJB, but crucially, there is no optimization involved.

This is a macro-labor model, where the amount of job openings depends on the value of a filled position ($V$), which is why the law of motion for $u$ (unemployment) depends on the equilibrium value of $V$.

How could I compute (numerically?) the value of $V$ that solves this ODE?

• Could you give further indications about $\lambda$ and other parameters ? Maybe it is better to see the problem by a dynamic optimization level with Hamiltonian. After, the value function could be numerically solved from it. Feb 1 '16 at 15:42