Seemingly unrelated error message: In an assignment A(I) = B, the number of elements in B and I must be the same.

Hello, I am trying to solve a ODE BVP using bvp4c. I cannot figure out why I am getting this error:

In an assignment  A(I) = B, the number of elements in B and I must be the same.

Error in HIVTBeqns (line 47)
ydot(1)=Lambda-mu.*X_SS-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_SS./N)-tau.*(X_ST+X_HT).*(X_SS./N);

Most of the following code consists of equations - the error message seems unrelated to what I have done!

%System of differential equations
function ydot = HIVTBeqns(t,y)

global Lambda mu mu_H mu_T mu_IRIS beta tau lambda_H lambda_H_1 lambda_H_2 lambda_T lambda_T_1 lambda_T_2 gamma psi alpha A B C D

X_SS = y(1);
X_ST = y(2);
X_HS = y(3); 
X_HT = y(4);
X_S_T_T_1 = y(5); 
X_S_T_T_2 = y(6);
X_S_T_T_3 = y(7);
X_H_T_T_1 = y(8);
X_H_T_T_2 = y(9); 
X_H_T_T_3 = y(10); 
X_H_T_T_T_1 = y(11); 
X_H_T_T_T_2 = y(12);
X_H_T_1 = y(13);
X_H_T_2 = y(14);
X_IRIS = y(15);
lambda_1 = y(16);
lambda_2 = y(17);
lambda_3 = y(18);
lambda_4 = y(19);
lambda_5 = y(20);
lambda_6 = y(21);
lambda_7 = y(22);
lambda_8 = y(23);
lambda_9 = y(24);
lambda_10 = y(25);
lambda_11 = y(26);
lambda_12 = y(27);
lambda_13 = y(28);
lambda_14 = y(29);
lambda_15 = y(30);

%%Hamiltonian and 'N'
N = X_SS + X_ST + X_HS + X_HT + X_S_T_T_1 + X_S_T_T_2 + X_S_T_T_3 + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3 + X_H_T_T_T_1 + X_H_T_T_T_2 + X_H_T_1 + X_H_T_2 + X_IRIS;

H = C.*(X_IRIS.*mu_IRIS + mu_H.*(X_H_T_T_1 + X_H_T_T_2)) - lambda_6.*(X_S_T_T_2.*(lambda_T_2 + mu) - X_S_T_T_1.*lambda_T_1 + (X_S_T_T_2.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N) - lambda_2.*(X_ST.*(mu + mu_T + lambda_T*psi) + (X_ST.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N - (X_SS.*tau.*(X_HT + X_ST))./N) - lambda_3.*(X_HS.*(lambda_H + mu + mu_H) - (X_SS.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N + (X_HS.*tau.*(X_HT + X_ST))./N) - lambda_7.*(X_S_T_T_3*mu - X_S_T_T_2*lambda_T_2 + (X_S_T_T_3*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N) + lambda_9.*(X_H_T_T_1*lambda_T_1 - X_H_T_T_2.*(lambda_H_2 + 2*lambda_T_2 + mu + mu_H) + (X_S_T_T_2.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N) - lambda_14.*(X_H_T_2.*(mu + (tau.*(X_HT + X_ST))./N) - X_HS.*lambda_H) + lambda_13.*(X_IRIS.*alpha + X_H_T_T_3.*lambda_H - X_H_T_1.*mu + lambda_T_2.*(X_H_T_T_2 + X_H_T_T_T_2)) - lambda_15.*(X_IRIS.*alpha - gamma.*(X_H_T_T_T_1 + (X_H_T_T_T_2.*lambda_T_1)/(lambda_T_1 + lambda_T_2)) + X_IRIS.*(mu + mu_IRIS)) + lambda_12.*(X_H_T_T_2.*lambda_H_2 + X_H_T_T_T_1.*lambda_T_1 - X_H_T_T_T_2.*(lambda_T_2 + mu) - (X_H_T_T_T_2.*gamma*lambda_T_1)/(lambda_T_1 + lambda_T_2)) + lambda_10.*(X_H_T_T_2.*lambda_T_2 - X_H_T_T_3.*(lambda_H + mu + mu_H) + (X_S_T_T_3.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N) + lambda_11.*(X_H_T_T_1.*lambda_H_1 - X_H_T_T_T_1.*(gamma + lambda_T_1 + mu) + (X_H_T_2.*tau.*(X_HT + X_ST))./N) - lambda_5.*(X_S_T_T_1.*(lambda_T_1 + mu) - X_ST.*lambda_T.*psi + (X_S_T_T_1.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N) + D.*(lambda_H_1^2 + lambda_H_2^2) + lambda_4.*((X_ST.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N - X_HT.*(lambda_T + mu + mu_H + mu_T) + (X_HS.*tau.*(X_HT + X_ST))./N) - lambda_1.*(X_SS*mu - Lambda + (X_SS*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N + (X_SS.*tau.*(X_HT + X_ST))./N) + lambda_8.*(X_HT.*lambda_T - X_H_T_T_1.*(lambda_H_1 + lambda_T_1 + mu + mu_H) + (X_S_T_T_1.*beta.*(X_HS + X_HT + X_H_T_T_1 + X_H_T_T_2 + X_H_T_T_3))./N) + B.*gamma.*(X_H_T_T_1 + (X_H_T_T_2.*lambda_T_1)./(lambda_T_1 + lambda_T_2)) + (A.*beta.*(X_H_T_T_1 + X_H_T_T_2).*(X_SS + X_ST + X_S_T_T_1 + X_S_T_T_2 + X_S_T_T_3))./N;

%%Optimality equations
lambda_H_1 = ((lambda_8-lambda_11).*X_H_T_T_1)./(2*D);
lambda_H_2 = ((lambda_9-lambda_12).*X_H_T_T_2)./(2*D);

%%Differential equations
ydot = zeros(30,1);
ydot(1)=Lambda-mu.*X_SS-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_SS./N)-tau.*(X_ST+X_HT).*(X_SS./N);
ydot(2)=tau.*(X_ST+X_HT).*(X_SS./N)-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_ST./N)-(mu+mu_T+psi.*lambda_T).*X_ST;
ydot(3)=beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_SS./N)-tau.*(X_ST+X_HT).*(X_HS./N)-(mu+lambda_H+mu_H).*X_HS;
ydot(4)=beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_ST./N)+tau.*(X_ST+X_HT).*(X_HS./N)-(mu+mu_H+mu_T+lambda_T).*X_HT;
ydot(5)=psi.*lambda_T.*X_ST-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_S_T_T_1./N)-(mu+lambda_T_1).*(X_S_T_T_1);
ydot(6)=lambda_T_1.*X_S_T_T_1-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_S_T_T_2./N)-(mu+lambda_T_2).*(X_S_T_T_2);
ydot(7)=lambda_T_2.*X_S_T_T_2-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_S_T_T_3./N)-(mu).*(X_S_T_T_3);
ydot(8)=lambda_T.*X_HT+beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_S_T_T_1./N)-(mu+mu_H+lambda_T_1+lambda_H_1).*X_H_T_T_1;
ydot(9)=lambda_T_1.*X_H_T_T_1+beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_S_T_T_2./N)-(lambda_H_2+mu+mu_H+lambda_T_2+lambda_T_2).*X_H_T_T_2;
ydot(10)=beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_S_T_T_3./N)+(lambda_T_2).*X_H_T_T_2-(mu+mu_H+lambda_H).*X_H_T_T_3;
ydot(11)=lambda_H_1.*X_H_T_T_1+tau.*(X_ST+X_HT).*X_H_T_2./N-(mu+lambda_T_1+gamma).*X_H_T_T_T_1;
ydot(12)=lambda_H_2.*X_H_T_T_2-(lambda_T_2+mu).*X_H_T_T_T_2+(lambda_T_1).*(X_H_T_T_T_1)-gamma.*(lambda_T_1./(lambda_T_1+lambda_T_2)).*(X_H_T_T_T_2);
ydot(13)=lambda_T_2.*(X_H_T_T_2+X_H_T_T_T_2)-mu.*X_H_T_1+alpha.*X_IRIS+(lambda_H).*X_H_T_T_3;
ydot(14)=lambda_H.*X_HS-(tau.*(X_ST+X_HT)./N+mu).*X_H_T_2;
ydot(15)=gamma.*(X_H_T_T_T_1+(lambda_T_1./(lambda_T_1+lambda_T_2)).*(X_H_T_T_T_2))-(mu+mu_IRIS).*X_IRIS-alpha.*X_IRIS;
ydot(16)=-diff(H,X_SS);
ydot(17)=-diff(H,X_ST);
ydot(18)=-diff(H,X_HS);
ydot(19)=-diff(H,X_HT);
ydot(20)=-diff(H,X_S_T_T_1);
ydot(21)=-diff(H,X_S_T_T_2);
ydot(22)=-diff(H,X_S_T_T_3);
ydot(23)=-diff(H,X_H_T_T_1);
ydot(24)=-diff(H,X_H_T_T_2);
ydot(25)=-diff(H,X_H_T_T_3);
ydot(26)=-diff(H,X_H_T_T_T_1);
ydot(27)=-diff(H,X_H_T_T_T_2);
ydot(28)=-diff(H,X_H_T_1);
ydot(29)=-diff(H,X_H_T_2);
ydot(30)=-diff(H,X_IRIS);
end

Here are the function files:

function res = HIVTBbcs(yb)

res = [yb(16)
       yb(17)
       yb(18)
       yb(19)
       yb(20)
       yb(21)
       yb(22)
       yb(23)
       yb(24)
       yb(25)
       yb(26)
       yb(27)
       yb(28)
       yb(29)
       yb(30)];
end

and

function solution
  solinit = bvpinit(linspace(0,1,30),[192941151 132478723 2422951 1817998 161501 161501 145865 7953 7953 117 3445 3445 3032 152450 82 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000]);
  sol = bvp4c(@HIVTBeqns1,@HIVTBbcs,solinit);
  end

Please let me know if I need me to fix anything in my presentation here.

I HAVE DONE

size (ydot(1)) and 

size (RHS), and am getting 

1 x 7 and 

1 x 1 respectively.

[Star Strider]

Everything on the RHS of ydot(1) (and the others) have to return scalars, the value of the particular derivative at that point. Something on the RHS is a vector, and that is likely causing the problem. You will have to go through the expression for ydot(1) (and the others as well), see what their components are, and be sure none are vectors.

(I’m not entering this as an Answer because it it’s only an observation.)

[AM]

Does the "diff" in

diff(H,X_SS)

etc. have anything to do with the problem? I have made sure everything is properly defined (.*, scalars not vectors, etc.)

[dpb]

Not likely, no, other than it implies X is a vector which may be a klew...see Answer.

[Star Strider]

Your ydot(1) (and the other individual elements of your ydot array) should each be (1x1). That ydot(1) is (1x7) won’t work. The problem is X_SS, X_HS, and the other vectors. You need to figure out what you want to do with the vectors, and how to choose a scalar from their values.

The diff statement is likely the problem, not in itself but because of its vector arguments. You need to figure out what you want with the diff statements as well, since you’ll get the same error with ydot(16), etc.

Also, if you want to take the numerical derivative, gradient is a better choice than diff.

Again, I’m not offering an Answer because I don’t have a solution to your problem.

Answer 1

[dpb]

In an assignment  A(I) = B, the number of elements in B and I must be the same.

Error in HIVTBeqns (line 47)

ydot(1)=Lambda-mu.*X_SS-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_SS./N)-tau.*(X_ST+X_HT).*(X_SS./N);

ydot(1) is a single array value/location whereas you've indicated above by the use of diff on X that it is a vector. That's also indicated by the use of the "dot" operators in the above that size will not return the values [1,1] in the operation

size(Lambda-mu.*X_SS-beta.*(X_HS+X_HT+X_H_T_T_1+X_H_T_T_2+X_H_T_T_3).*(X_SS./N)-tau.*(X_ST+X_HT).*(X_SS./N))

That'll be the problem, specifically; now what the actual needed size will be is up to your problem statement and I've not tried to wade thru the amount of code.

[AM]

Ok. I have rewritten all the "diff" as actual equations, with all of the * as .* and \ as .\ Also, length(H)=length(Lambda). But I am still getting the same error message!

[dpb]

Yes, to assign to ydot(1) you'll have to have the RHS be one, single value, no vectors in sight.

The question IA and I can't answer is how to recast the problem statement into the proper form, specifically, only what the syntax error is and means.