I’m just messaging to ask if there is any plans to implement an adaptive mesh in the ODE-BVP solvers in BoundaryValueDiffEq.jl? I’m working on a problem where this is necessary, and for which the adaptive mesh in bvp4c in Matlab works quite well. Even the ability to specify a non-uniform mesh (rather than only specifying the dt value for a uniform mesh) would be very helpful.
Thanks
It’s planned… we just never got back to it. To teach @yingboma how adaptivity and interpolations of ODEs work, I had him do some ODE stuff and then he started optimizing our stiff ODE solvers to benchmark super well haha. We should return to that soon since it’s really not difficult and it’s just waiting to be done… (TBH the main reason is probably because no one on the dev team has very many BVPs to solve, so if you give a good reason we might whip adaptivity together 
)
We have tons of them (almost 100?). This is specifically asking for an adaptive MIRK method for BVPs.
Isn’t that what approxfun is good at?
Haha, it’s good to hear the backstory. DifferentialEquations.jl is incredible - thanks for all your work. I thought of implementing adaptivity myself - as you say, conceptually it’s straightforward enough. However, a brief look through the source code for DifferentialEquations.jl has reminded me I might be more of a Julia user than a developer 
You asked for a planned use-case, so…
In chemical engineering, BVPs are arguably more common than IVPs. This is because we’re often interested in modelling spatially non-uniform physical systems, and in such cases you tend to get BCs in multiple locations simultaneously. For example, when designing distillation columns or absorbers, we’ll often specify some concentrations at the top and some at the bottom, and solve for the concentrations everywhere in between. Similarly, when dealing with reactive absorption of a gas into a liquid, we’ll have some boundary conditions at the surface of the liquid, and some from the bulk of the liquid. Secretly, I’d like to combine both into one complete model - right now such multi-scale models (simultaneously modelling a whole column and also of the local gas-liquid reactive dynamics) are only available in expensive, commercial solvers, where special care is taken to discretise everything just right and form carefully-structured matrices which are solved via Fortran. I suspect a julia implementation built off DifferentialEquations.jl could be conceptually much simpler and fast enough to be of practical use.
If it’s a smooth BVP, then yes ApproxFun.jl is a good option, and probably the best option.
That’s an interesting thought. Can approxfun handle stiff systems? I’ve found that Matlab’s bvp4c is particularly capable in that regard.
Yeah, and that’s the part we cared to do. It has special forms for two-point BVPs so that it can utilize BandedMatrix types and then exploit the banded QR algorithm in that case.
Anything that isn’t shooting should do well there, since BVPs are inherently implicit. If your problem is ill-conditioned and smooth, spectral discretizations probably do better than MIRK methods.
Thanks for your help - I will give ApproxFun,jl a try.
Just to follow up on this - do you know if it is possible to define domains with breakpoints in ApproxFun.jl, as discussed here for Chebfun:
https://www.chebfun.org/examples/ode-linear/Breakpoints.html
This would probably be ‘close enough’ to adaptive meshing for boundary layer problems. I’ve played around with the various ‘PiecewiseSegment’ options, but they don’t seem to play nice with the newton solver in ApproxFun.jl.
Thanks
In theory it should work, though it’s combining two “proof of concept” pieces, so will require a number of bug fixes. Unfortunately I’m in the middle of completely revamping the underpinnings of ApproxFun so that it scales better (including to fully support hp methods) so I won’t have time to look at this, though any PR will be greatly appreciated.
Which two proof of concept pieces would need to be combined? I might have a look - I understand the theory, but don’t really understand Julia
Hello, I also would like to use a specified mesh rather than just defining a dt. Did this update already happen? If so, could you point me to documentation on it? Thanks!
I’m not sure it takes a vector dt at this time.
There are many different BVP solvers in BifurcationKit.jl. However, they are tailored for periodic orbits. Changing the boundary condition would not take too much effort. For example, you would need to change this line for the collocation method.
You can specify the time mesh for all of them.
You can automatically adapt the time mesh for the orthogonal collocation algo during continuation. I would say, with a newton callback, you could perhaps adapt during newton iterations.
However, it is not ready for use as is for your use case.