Solving fully implicit linear boundary value problem with constants that change along the domain

I have a system of implicit linear steady state ODEs with 8 unknown variables () and 8 equations to be solved along r. All the sets of equations contain first order derivatives only and known coefficient terms (). The coefficient terms also change with r, however, the values of each of these terms at every point in r are known and already stored as MATLAB vectors. The equation sets are of the form:

.

.

.

I have the following boundary conditions

I attempted to use the mass matrix method https://www.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html to reduce it until I got a function handle that could be passed to an ODE solver. This cannot work however, because the method seems suited to IVPs or BVPs where all the boundary conditions are defined at the beginning of the domain only. Also, supposing I figure out a way to define all my boundary conditions at the beginning of the domain, I can't still figure out a way to use decic to get consistent set of initial conditions because my coefficient terms are vectors and not scalar variables, since they also change with r and are not constants..

Error using symengine
Array sizes must match.

Error in sym/privBinaryOp (line 1032)
            Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in  +  (line 7)
X = privBinaryOp(A, B, 'symobj::zipWithImplicitExpansion', '_plus');

Error in
symengine>@(r,in2,in3,param14,param15,param16,param17,param18,param19,param20,param21,param22,param23,param24,param25,param26)[in2(1,:).*(param22+param26./r)+in2(3,:).*(param20+param16./r)+in3(3,:).*param16+in3(1,:).*param26-(in2(4,:).*param17.*param25)./r-(in2(6,:).*param25.*param26)./r;in2(2,:).*(param22+param26./r)+in2(4,:).*(param20+param16./r)+in3(4,:).*param16+in3(2,:).*param26+(in2(3,:).*param17.*param25)./r+(in2(5,:).*param25.*param26)./r;in3(7,:)-in2(3,:).*(-param16.^2./r+param17.^2./r-param16.*param20.*2.0+param17.*param23.*2.0+param23.^2.*r)-in2(5,:).*(param23.*param26.*2.0+(param17.*param26.*2.0)./r)+in3(3,:).*param16.^2+in2(1,:).*(param16.*param22.*2.0+param20.*param26.*2.0+(param16.*param26.*2.0)./r)+in3(1,:).*param16.*param26.*2.0-(in2(2,:).*param17.*param25.*param26)./r-(in2(4,:).*param16.*param17.*param25)./r-(in2(6,:).*param16.*param25.*param26)./r;in3(8,:)-in2(4,:).*(-param16.^2./r+param17.^2./r-param16.*param20.*2.0+param17.*param23.*2.0+param23.^2.*r)-in2(6,:).*(param23.*param26.*2.0+(param17.*param26.*2.0)./r)+in3(4,:).*param16.^2+in2(2,:).*(param16.*param22.*2.0+param20.*param26.*2.0+(param16.*param26.*2.0)./r)+in3(2,:).*param16.*param26.*2.0+(in2(1,:).*param17.*param25.*param26)./r+(in2(3,:).*param16.*param17.*param25)./r+(in2(5,:).*param16.*param25.*param26)./r;-in2(8,:).*param25+in2(1,:).*(param17.*param26.*2.0+param17.*param22.*r+param21.*param26.*r+param23.*param26.*r.*2.0)+in2(3,:).*(param16.*param17.*2.0+param16.*param20.*r+param17.*param20.*r+param17.*param23.*r.*2.0)+in2(5,:).*(param16.*param26.*2.0+param16.*param22.*r+param20.*param26.*r)-in2(4,:).*param17.^2.*param25-in2(6,:).*param17.*param25.*param26.*2.0+in3(3,:).*param16.*param17.*r+in3(1,:).*param17.*param26.*r+in3(5,:).*param16.*param26.*r;-in2(7,:).*param25+in2(2,:).*(param17.*param26.*2.0+param17.*param22.*r+param21.*param26.*r+param23.*param26.*r.*2.0)+in2(4,:).*(param16.*param17.*2.0+param16.*param20.*r+param17.*param20.*r+param17.*param23.*r.*2.0)+in2(6,:).*(param16.*param26.*2.0+param16.*param22.*r+param20.*param26.*r)+in2(3,:).*param17.^2.*param25+in2(5,:).*param17.*param25.*param26.*2.0+in3(4,:).*param16.*param17.*r+in3(2,:).*param17.*param26.*r+in3(6,:).*param16.*param26.*r;in2(3,:).*((param20./2.0+param16./(r.*2.0)).*(-param23.^2.*r.^2+param16.^2+param17.^2)+param16.*(param16.*param20+param17.*param21-param23.^2.*r))+in3(1,:).*(param14+param15+param16.^2.*param26)+in2(5,:).*(param16.*param17.*param22+param16.*param21.*param26+param17.*param20.*param26+(param16.*param17.*param26)./r)+in2(1,:).*(param18+param19+param16.^2.*param22+param14./r+param15./r+(param16.^2.*param26)./r+param16.*param20.*param26.*2.0)+(in3(3,:).*param16.*(-param23.^2.*r.^2+param16.^2+param17.^2))./2.0-in2(6,:).*param25.*(param14./r+param15./r+(param17.^2.*param26)./r)+in3(5,:).*param16.*param17.*param26+(in2(7,:).*param24.*(param20+param16./r))./(param24-1.0)+(in3(7,:).*param16.*param24)./(param24-1.0)-(in2(4,:).*param17.*param25.*(-param23.^2.*r.^2+param16.^2+param17.^2))./(r.*2.0)-(in2(8,:).*param17.*param24.*param25)./(r.*(param24-1.0))-(in2(2,:).*param16.*param17.*param25.*param26)./r;in2(3,:).*((param20./2.0+param16./(r.*2.0)).*(-param23.^2.*r.^2+param16.^2+param17.^2)+param16.*(param16.*param20+param17.*param21-param23.^2.*r))+in3(2,:).*(param14+param15+param16.^2.*param26)+in2(6,:).*(param16.*param17.*param22+param16.*param21.*param26+param17.*param20.*param26+(param16.*param17.*param26)./r)+in2(2,:).*(param18+param19+param16.^2.*param22+param14./r+param15./r+(param16.^2.*param26)./r+param16.*param20.*param26.*2.0)+(in3(4,:).*param16.*(-param23.^2.*r.^2+param16.^2+param17.^2))./2.0+in2(5,:).*param25.*(param14./r+param15./r+(param17.^2.*param26)./r)+in3(6,:).*param16.*param17.*param26+(in2(8,:).*param24.*(param20+param16./r))./(param24-1.0)+(in3(8,:).*param16.*param24)./(param24-1.0)+(in2(3,:).*param17.*param25.*(-param23.^2.*r.^2+param16.^2+param17.^2))./(r.*2.0)+(in2(7,:).*param17.*param24.*param25)./(r.*(param24-1.0))+(in2(1,:).*param16.*param17.*param25.*param26)./r]

Error in Equations (line 100)
F = @(r,Y,YP) f(r,Y,YP,E,P_bar,Vr_bar,Vt_bar,dE,dP_bar,dVr_bar,dVt_bar,drho_bar,ff,...

Error in decic (line 66)
res = feval(odefun,t0,y0,yp0,varargin{:});

Error in Equations (line 106)
[y0,yp0] = decic(F,0,y0est,[],yp0est,[],opt)

1. How can I use decic to get a set of consistent initial conditions, despite having variable coefficients that depend on r (if I decide to go ahead with this idea and set all my BCs to be defined at the beginning of the domain)?
2. Is there a better way to solve an implicit set of linear steady state ODEs with all kinds of boundary conditions (Dirichlet, Neumann) in MATLAB?