Ok, so looks like what you are after is solving a DAE that has BCs at more than one boundary.
This system seems simple enough to use a shooting strategy, I think. You have 3 different boundary locations, so you can treat the problem as if it were an IVP leaving 2 guesses for 2 out of 3 of the Y's at x=0 as guesses (not sure which ones, you are mixing a's, b's, and 1,2,3's), use either of the odeXX (with appropriate handling of the algebraic equation), and you will get some values for the Y's at the 2 other boundaries in the course of integration. Their deviations from your actual BCs will become residuals for a 2-variable root finding.
Best Answer