MATLAB: Solving BVP for DAE in Matlab

bvpbvp4cbvp5cdaeodeode15i

How do I solve boundary value problems (BVP) for differential algebraic equation systems (DAE) in Matlab?
bvp4c and bvp5c can only handle differential equations. ode15i can handle DAE, but only initial value problems (IVP).
Transformation of the DAE into a differential equations system is too costly or impossible.
Also, I do not want to write my own BVP numerical routine that uses matlabs differential equation solvers (we use Matlab because when we do not want to write our own routines, right?).
Currently I am using bvp4c/bvp5c and the function containing the differential equations calls another functions that solved the algebraic (non-differential) subsystem before calculating the rhs of the differential equations. However, this seems highly unstable and inefficient to me.
Any recommendations or experience on solving BVP for DAE in Matlab is highly welcome.
Best,
Benjamin