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.

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])

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.

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.

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) 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?

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

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.