Did matlab really fail to solve this system of 4 DAEs or did I make a mistake somewhere?

Hi everyone, I'm trying to analytically solve a system of 4 DAEs, however all outputs are blank. Please see the code below:

function WaterLiquidFilmProfilesAnalytical

syms a b c e f t A B C E F G H K

dEq1 = 'A * D2a + B * D2b + C * D2c = 0';

dEq2 = 'E * D2e - B * D2b - 2 * C * D2c - F * D2f = 0';

      Eq1 = ('G * a = e * b'); 

              Eq2 = str2sym(Eq1);

                  [dEq3, initEq3] = TurnEqIntoDEq(Eq2, [G a e b], t, 0);

dEq3_char = SymArray2CharCell(dEq3);
                            initEq3_char = SymArray2CharCell(initEq3);

        Eq4 = ('H * b = e * c'); 

              Eq5 = str2sym(Eq4);

                  [dEq4, initEq4] = TurnEqIntoDEq(Eq5, [H b e c], t, 0);

dEq4_char = SymArray2CharCell(dEq4);
                            initEq4_char = SymArray2CharCell(initEq4);

              Eq6 = ('K = e * f'); 

                      Eq7 = str2sym(Eq6);

                  [dEq5, initEq5] = TurnEqIntoDEq(Eq7, [K e f], t, 0);

dEq5_char = SymArray2CharCell(dEq5);
                            initEq5_char = SymArray2CharCell(initEq5);

%[sol_dEq1, sol_dEq2, sol_dEq3, sol_dEq4, sol_dEq5] = dsolve(dEq1, dEq2, dEq3_char{:},'a(0)=0','e(0)=0','b(0)=0', initEq3_char{:}, 't', dEq4_char{:},'b(0)=0','e(0)=0','c(0)=0', initEq4_char{:}, 't', dEq5_char{:},'d(0)=0','f(0)=0', initEq5_char{:}, 't')

%[sol_dEq1, sol_dEq2, sol_dEq3, sol_dEq4, sol_dEq5] = dsolve(dEq1, dEq2, dEq3_char{:},  dEq4_char{:},  dEq5_char{:},'a(0)=0','b(0)=0','c(0)=0','e(0)=0','f(0)=0',initEq3_char{:},initEq4_char{:}, initEq5_char{:},'t')

[a, b, c, d, e] = dsolve(dEq1, dEq2, dEq3_char{:},  dEq4_char{:},  dEq5_char{:},'a(0)=0','b(0)=0','c(0)=0','e(0)=0','f(0)=0',initEq3_char{:},initEq4_char{:}, initEq5_char{:},'t')
end

function [D_Eq, initEq] = TurnEqIntoDEq(eq, depVars, indepVar, initialVal)
% eq = equations
% depVars = dependent variables
% indepVar = independent variable
% initialVal = initial value of indepVar

depVarsLong = sym(zeros(size(depVars)));
for k = 1:numel(depVars)
  % Making the variables functions
  % eg. a becomes a(t)
  % This is so that diff(a, t) does not become 0
  depVarsLong(k) = str2sym([char(depVars(k)) '(' char(indepVar) ')']);
end

% Next making the equation in terms of these functions
eqLong = subs(eq, depVars, depVarsLong);

% Now find the ODE corresponding to the equation
D_EqLong = diff(eqLong, indepVar);

% Now replace all the long terms like 'diff(a(t), t)'
% with short terms like 'Da'

D_depVarsShort = sym(zeros(size(depVars)));
for k = 1:numel(depVars)
  D_depVarsShort(k) = str2sym(['D' char(depVars(k))]);
end
% Next make the long names like 'diff(a(t), t)'
D_depVarsLong = diff(depVarsLong, indepVar);
% Finally replace
D_Eq = subs(D_EqLong, D_depVarsLong, D_depVarsShort);

% Finally determine the equation
% governing the initial values
initEq = subs(eqLong, indepVar, initialVal);
end

function cc = SymArray2CharCell(sa)
cc = cell(size(sa));
for k = 1:numel(sa)
  cc{k} = char(sa(k));
end

end

Did matlab really fail or did I make a mistake?

I also tried to solve these DAEs numerically,using:

syms a(x) b(x) c(x) d(x) e(x) A B C D E FF G H 

eqn1 = A * diff(a(x),x,2) + B * diff(b(x),x,2) + C * diff(c(x),x,2) == 0;

eqn2 = D * diff(d(x),x,2) - B * diff(b(x),x,2) - 2 * C * diff(c(x),x,2) - E * diff(e(x),x,2) == 0;

eqn3 = FF == (d(x) * b(x)) / a(x);

eqn4 = G == (d(x) * c(x)) / b(x);

eqn5 = H == d(x) * e(x);

eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];

vars = [a(x); b(x); c(x); d(x); e(x)];

origVars = length(vars);

M = incidenceMatrix(eqns, vars);

[eqns, vars] = reduceDifferentialOrder(eqns, vars);

isLowIndexDAE(eqns,vars);

f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H);

A = 2.89e-9;   
B = 2.35e-9;  
C = 1.69e-9;  
D = 1.6e-8;    
E = 9.25e-9; 
FF = 6.24;      
G = 5.68e-5;  
H = 5.3e-8; 

F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H);

vars;

y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];
yp0est = zeros(5,1);
opt = odeset('RelTol', 10^(-3), 'AbsTol' , 10^(-3));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 7.21e-05], y0, yp0, opt);
for k = 1:origVars
  S{k} = char(vars(k));
end

plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
legend('SO_2','HSO_{3}^{-}', 'SO_{3}^{2-}', 'H^+','OH^-','location','Best')

However, I get the error message:

Index in position 1 exceeds array bounds (must not exceed 5).

Error in
symengine>@(x,in2,in3,param1,param2,param3,param4,param5,param6,param7,param8)[in3(6,:).*param1+in3(7,:).*param2+in3(8,:).*param3;-in3(7,:).*param2-in3(8,:).*param3.*2.0+in3(9,:).*param4-in3(10,:).*param5;param6-(in2(2,:).*in2(4,:))./in2(1,:);param7-(in2(3,:).*in2(4,:))./in2(2,:);param8-in2(4,:).*in2(5,:);in2(6,:)-in3(1,:);in2(7,:)-in3(2,:);in2(8,:)-in3(3,:);in2(9,:)-in3(4,:);in2(10,:)-in3(5,:)]

Error in WaterProfiles>@(t,Y,YP)f(t,Y,YP,A,B,C,D,E,FF,G,H)

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

Error in WaterProfiles (line 45)
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

Debugging points me to line 45:

[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

And in my understanding, the index of y0est does not exceed 5:

y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];

[Stephan]

This appears to be the same question like here:

https://de.mathworks.com/matlabcentral/answers/417194-what-is-the-meaning-of-index-in-position-1-exceeds-array-bounds-must-not-exceed-5

[Dursman Mchabe]

Hi Stephan. The question on the link does not have the analytical (dsolve) code.

[Stephan]

Is this question anwsered with the answer of the other one? The cause of the error message should now have been clarified. Or is there more to do here?

[Dursman Mchabe]

Hi Stephen, thank you so much for providing the numerical solution. I truly appreciate. I have assumed that the analytical solution is impossible. Is my assumption correct?

[Stephan]

i dont know - i did not try. I just asked because the error message is the same then in your other question.

[Dursman Mchabe]

My wish was to get an analytical solution. For consistency with other sections of my work. Because, if I use a numerical solution, I will be forced to solve a three-point bvp problem. Which is very very difficult.

[Stephan]

Ok. Then i guess the question is still open. Maybe someone can help.