# Singular BVP with integrable boundaries

**URL:** <https://discourse.julialang.org/t/singular-bvp-with-integrable-boundaries/135224>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, differentialequation, ordinarydiffeq, bvp\
**Created:** [January 23, 2026, 2:36pm UTC](https://discourse.julialang.org/t/singular-bvp-with-integrable-boundaries/135224 "2026-01-23T14:36:53Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![NoFishLikeIan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nofishlikeian/32/20917_2.png) [@NoFishLikeIan](https://discourse.julialang.org/u/NoFishLikeIan)\
**Post date:** [January 23, 2026, 2:36pm UTC](https://discourse.julialang.org/t/singular-bvp-with-integrable-boundaries/135224/1 "2026-01-23T14:36:53Z")

</div>

I’m solving a singular first-order, two dimensional ODE arising in reputation games ([Faingold & Sannikov, 2011](https://doi.org/10.3982/ECTA7377)).

On \phi \in [0, 1], I have

u' = \frac{v}{\phi(1-\phi)} \\ v' = \frac{1}{\phi(1-\phi)} \left(v + 2\rho \frac{(u - w(\phi, v/\rho))}{\gamma(v / \rho)^2} \right)

with boundary conditions v(0) = v(1) = 0, u(0) = w(1, 0) and u(1) = w(0, 0). Here \gamma(\cdot) \> 0. Hence, at the boundaries, I have integrable singularities.

My current approach has been to try and solve the problem numerically using a `TwoPointBVProblem` in the interior \phi \in [\varepsilon, 1 - \varepsilon] for some small \varepsilon, using the `MIRK4` method. But unfortunately, I get instability.

Can anyone suggest an algorithm or a change of variable to solve this?

I am adding a simplified **MWE** of the problem below, to illustrate my current approach.

```julia
using DifferentialEquations

function w(ϕ, z)
    0.1 * ϕ + 0.05 * z
end

function γ(z)
    1.0 + z
end

function F!(dx, x, ρ, φ)
    u, v = x
    
    dx[1] = v / (φ * (1 - φ))
    dx[2] = (v + 2ρ * (u - w(φ, v / ρ)) / γ(v / ρ)^2) / (φ * (1 - φ))
end

function leftbc!(res, x₀, p)
    u₀, v₀ = x₀
    res[1] = u₀ - w(0., 0.)
    res[2] = v₀
end

function rightbc!(res, x₁, p)
    u₁, v₁ = x₁
    res[1] = u₁ - w(1., 0.)
    res[2] = v₁
end

ρ = 0.05
x₀ = [w(0., 0.), 0.]

ε = 1e-2
bvp = TwoPointBVProblem(F!, (leftbc!, rightbc!), x₀, (ε, 1 - ε), ρ; bcresid_prototype = (zeros(2), zeros(2)))
sol = solve(bvp, MIRK4(), dt = 0.01)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 24, 2026, 7:01pm UTC](https://discourse.julialang.org/t/singular-bvp-with-integrable-boundaries/135224/2 "2026-01-24T19:01:38Z")

</div>

Open an issue, we don’t handle this case so far.

---

<div class="post-metadata">

**Author:** ![NoFishLikeIan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nofishlikeian/32/20917_2.png) [@NoFishLikeIan](https://discourse.julialang.org/u/NoFishLikeIan)\
**Post date:** [January 26, 2026, 7:13am UTC](https://discourse.julialang.org/t/singular-bvp-with-integrable-boundaries/135224/3 "2026-01-26T07:13:46Z")

</div>

Thanks, will do. Do you have a hunch/reference on how to solve these? If I can implement it myself, I can contribute to the issue.