# Solving a second order coupled BVP with boundaries at infinity using DifferentialEquations.jl

**URL:** https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168
**Category:** New to Julia
**Tags:** question, diffeq
**Created:** [April 24, 2020, 10:40pm UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168 "2020-04-24T22:40:46Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![JonathanLoz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonathanloz/32/14311_2.png) [@JonathanLoz](https://discourse.julialang.org/u/JonathanLoz)
#### Post date: [April 24, 2020, 10:40pm UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168/1 "2020-04-24T22:40:46Z")

</div>

I’m trying to solve the following coupled differential equations in Julia by using DifferentialEquations.jl and BoundaryValueDiffEq.jl

\begin{align} \frac{d^2f}{du^2}+\frac{1}{u}\frac{df}{du}-\frac{n^2}{u^2}(1-a)^2f-\frac{\lambda}{e^2}(f-1)(f-\tilde{V}\_1)(f-\tilde{V}\_m)&=0\\ \frac{d^2a}{du^2}-\frac{1}{u}\frac{da}{du}+2f^2(1-a)&=0 \end{align}

With boundary conditions a[0]=f[0]=0, a[\infty]=f[\infty]=1.

Indeed, I already got some solutions using the GeneralMIRK4() solver.The problem arise when I increase the interval of integration and the compilation takes for ever to end. This might be happening because I’m not implementing the best method to solve the equations. Any advice in how to speed up my code or how to implement another method to solve this system of differential equations is appreciated.

```julia
function vortexgg!(dx,x,p,t)
    dx[1] = x[2]
    dx[2] = -x[2]/t+x[1]*(1.0-x[3])^2*(1/t^2)+
    (14.0)*(x[1]-1)*(x[1]-(0.3))*(x[1]-(0.6))

    dx[3] = x[4]
    dx[4] = x[4]/t-(2.0*x[1]^2)*(1.0-x[3])
 end

 function bc1!(residual, u, p, t)
    residual[1] = u[1][1] - 0.0
    residual[2] = u[end][1] - 1.0
    residual[3] = u[1][3] - 0.0
    residual[4] = u[end][3] - 1.0

 end
 bvp1 = BVProblem(vortexgg!, bc1!, [0.00,0.0,0.0,0.00], tspan)
  st = solve(bvp1, GeneralMIRK4(),dt=0.1,alg_hints = [:stiff])

```

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 25, 2020, 12:54am UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168/2 "2020-04-25T00:54:30Z")

</div>

Try doing a change of variables like u = t/(1-t) that maps [0,\infty) to [0,1).

---

<div class="post-metadata">

### Author: ![JonathanLoz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonathanloz/32/14311_2.png) [@JonathanLoz](https://discourse.julialang.org/u/JonathanLoz)
#### Post date: [April 25, 2020, 1:09am UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168/3 "2020-04-25T01:09:46Z")

</div>

Thanks, I already tried doing that but somehow it makes the system more instable, and the functions don’t converge to the appropiate solution. I’ll give it another try, I’m guessing that this could speed up the code.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 25, 2020, 5:40am UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168/4 "2020-04-25T05:40:37Z")

</div>

I am not sure what BoundaryValueDiffEq.jl uses, but I would consider simply using collocation on either [0,1] or [0, \infty] (which is a bit more tricky), see _Boyd, J. P. (2001). Chebyshev and fourier spectral methods._

---

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 25, 2020, 7:47am UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168/5 "2020-04-25T07:47:02Z")

</div>

> [@JonathanLoz](#):
>
> GeneralMIRK4

I am not sure this is a good choice on large intervals.

Do you know where `(f,a)` converge when `u` goes to infinity?

In any cases, I would give a go to `Shooting` (see [docs](https://docs.sciml.ai/stable/solvers/bvp_solve/index.html)) with appropriate ode stepper.

Finally, the factor `1/u` is not very nice from an ODE point of vue. Keep in mind that the solvers implement existence results from theorem. I see that your vector field is singular, a well chosen change of time variable could make it more well posed?

---

<div class="post-metadata">

### Author: ![bilderbuchi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bilderbuchi/32/13562_2.png) [@bilderbuchi](https://discourse.julialang.org/u/bilderbuchi)
#### Post date: [May 14, 2020, 7:32pm UTC](https://discourse.julialang.org/t/solving-a-second-order-coupled-bvp-with-boundaries-at-infinity-using-differentialequations-jl/38168/6 "2020-05-14T19:32:16Z")

</div>

So, this topic kinda piqued my interest and kept simmering in the back of my head. 😁 I finally sat down to work on an _asymptotic expansion_ of the problem at the boundaries.  
The behaviour of the functions, especially as u\rightarrow \infty in the present problem, can either guide you to an appropriate coordinate transformation (e.g. to make the functions decay nicely in the new coordinates), or you can use them to shift the right boundary to a location before \infty if you get a singularity there. I used both of these techniques to solve a BVP during my PhD work.

It’s been quite a while since I last did that, so I am pretty rusty. I started at u\rightarrow 0 first, and found

f\_L(u\rightarrow 0) \sim u\cdot f\_1 - u^2 \cdot f\_2 + \ldots\\ a\_L(u\rightarrow0) \sim u^2 \cdot a\_1 - u^4 \cdot a\_2 + \ldots

Is this useful to you? Values for the coefficients f\_1 and a\_1 are unknown, but that’s not the primary result of the analysis anyway. Interestingly, I have expressions for f\_2 and a\_2 (with some confidence).

As expected, the right boundary expansion turned out to be trickier, so I don’t have a result to show, yet – the effort passed my threshold for “nice evening fun”, the first part was ~8 pages of paper already.
