Index in position 2 exceeds array bounds

Could someone please tell me what should I put for D2q in odefcn?

syms p Dp D2p q Dq D2q x  
r=0.06;
s=0.3;
m=0.01;
f=0.075;                    
vmax=5;
smax=3;
xf=90;
w=0.15;
l=0.01;
a0=0.001;
h=0.2;
alpha=1;
pi=0.5;
mk=5;
mum=0.6;
t=0.1;
tp=0.05;
k=0.3;
a=((q+x*Dq)*t+(2*Dq*x+(x^2)*D2q)*s*tp)/(k-(2*Dq*x+(x^2)*D2q)*(tp^2));
m1=(Dp/p)*x*(m+a*t-mum)+0.5*(x^2)*((s+a*tp)^2)*(D2p/p);
t11=((w+m-r-f+mum)*(sqrt(vmax*alpha*xf)))/((pi^2)*sqrt(mk));
ode = p-(((x*alpha)*(w+m1+mum-r-f)^2)/((pi^4)));
D2ynum = solve(ode==0,D2p);
D2ynum = D2ynum(2);
f1 = matlabFunction(D2ynum,"Vars",{x, [p Dp q Dq D2q]})

f1 = function_handle with value:
    @(x,in2)in2(:,1).*1.0./x.^2.*1.0./((in2(:,3)./1.0e+1+in2(:,4).*x.*(1.3e+1./1.0e+2)+in2(:,5).*x.^2.*(3.0./2.0e+2))./((in2(:,4).*x)./1.0e+1+(in2(:,5).*x.^2)./2.0e+1-6.0)-3.0./1.0e+1).^2.*((sqrt(in2(:,1)).*1.0./sqrt(x))./4.0-(in2(:,2).*x.*((in2(:,3)./1.0e+1+in2(:,4).*x.*(1.3e+1./1.0e+2)+in2(:,5).*x.^2.*(3.0./2.0e+2))./((in2(:,4).*x)./2.0e+1+(in2(:,5).*x.^2)./4.0e+1-3.0)+5.9e+1./1.0e+2))./in2(:,1)+1.23e+2./2.0e+2).*-2.0

ode1=q-((((p*f/x)-mum*x*Dq+(m+a*t)*Dq*x+0.5*((s+a*tp)^2)*(2*x*Dq+(x^2)*D2q)-0.5*k*(a^2))/(r-(m+a*t))));
D2ynum1 = solve(ode1==0,D2q);

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

D2ynum1 = D2ynum1(1);
f2 = matlabFunction(D2ynum1,"Vars",{x, [p Dp q Dq]})

f2 = function_handle with value:
    @(x,in2)-(1.0./x.^2.*(in2(:,1).*1.8e+3-in2(:,3).*x.*1.2e+3-in2(:,4).*x.^2.*1.2e+4+in2(:,3).^2.*x.*4.0e+2+in2(:,4).^2.*x.^3.*8.76e+2-in2(:,4).*in2(:,1).*x.*3.0e+1+in2(:,4).*in2(:,3).*x.^2.*1.06e+3))./(in2(:,1).*-1.5e+1+x.*1.08e+3+in2(:,3).*x.*1.3e+2+in2(:,4).*x.^2.*2.38e+2)

odefcn = @(x,y)[y(2);f1(x, [y(1) , y(2), y(3), y(4)]);y(4);f2(x, [y(1), y(2), y(3) ,y(4)])];
bcfcn = @(ya,yb)[ya(1);yb(1)-t11;ya(4);yb(3)-((f*t11-0.5*k*(a^2)/((r-(m+t*a)))))];
xmesh = linspace(0.001,xf,10);
solinit = bvpinit(xmesh, [0.001 0.001 0.001 0.001]);
sol = bvp4c(odefcn,bcfcn,solinit);

Index in position 2 exceeds array bounds. Index must not exceed 4.

Error in symengine>@(x,in2)in2(:,1).*1.0./x.^2.*1.0./((in2(:,3)./1.0e+1+in2(:,4).*x.*(1.3e+1./1.0e+2)+in2(:,5).*x.^2.*(3.0./2.0e+2))./((in2(:,4).*x)./1.0e+1+(in2(:,5).*x.^2)./2.0e+1-6.0)-3.0./1.0e+1).^2.*((sqrt(in2(:,1)).*1.0./sqrt(x))./4.0-(in2(:,2).*x.*((in2(:,3)./1.0e+1+in2(:,4).*x.*(1.3e+1./1.0e+2)+in2(:,5).*x.^2.*(3.0./2.0e+2))./((in2(:,4).*x)./2.0e+1+(in2(:,5).*x.^2)./4.0e+1-3.0)+5.9e+1./1.0e+2))./in2(:,1)+1.23e+2./2.0e+2).*-2.0

Error in solution>@(x,y)[y(2);f1(x,[y(1),y(2),y(3),y(4)]);y(4);f2(x,[y(1),y(2),y(3),y(4)])] (line 31)
odefcn = @(x,y)[y(2);f1(x, [y(1) , y(2), y(3), y(4)]);y(4);f2(x, [y(1), y(2), y(3) ,y(4)])];

Error in bvparguments (line 105)
testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 122)
bvparguments(solver_name,ode,bc,solinit,options,varargin);

Answer 1

[Torsten]

I already showed you how to solve "ode" as well as "ode1" simultaneously for D2p and D2q. This is necessary because both depend on D2q.

Answer 2

[Walter Roberson]

f1 = matlabFunction(D2ynum,"Vars",{x, [p Dp q Dq D2q]})

5 variables bunched together, so f1 is going to expect to be able to index up to 5 for the second parameter.

odefcn = @(x,y)[y(2);f1(x, [y(1) , y(2), y(3), y(4)]);y(4);f2(x, [y(1), y(2), y(3) ,y(4)])];

but there you are only passing 4 values to the second parameter of f1.

[Della]

Thanks @Torsten and @Walter Roberson.