Subscript indices must either be real positive integers or logicals.

my main program is

Vd=40; %volume of drum m^3
Vr=37; %volume of riser m^3
Vdc=11; %volume of downcomer m^3
Vt=Vd+Vr+Vdc;%total volume
Ad=20; %drum area at normal operating level m^2
md=95000;%mass of drum in kg
mt=300000; %total metal mass kg
mr=160000; %total riser mass kg
k=25;%friction coefficent in downcomer riser loop
B=0.3;%parameter
Td=12; % residence time of steam min drum s
Cp=0.42;%specific heat of steel Kj/kgdegC
n=100;
t=zeros(1,n);
t(1,1)=0;
qf=50; %feed flow rate kg/s(to be specified by user)
p=zeros(1,n);
p(1,1)=8.5; %steam pressue MPa(to be specified by user)
Vwt=zeros(1,n);
Vwt(1,1)=57.2;%initial total volume of water in m^3
Vst=zeros(1,n);
Vst(1,1)=Vt-Vwt(1,1);
pf=9.142235;%pressure at which feed is entering Mpa(to be specified by user)
hf=1026.33;%-0.0003*pf^6+0.019*pf^5-0.47*pf^4+6.1*pf^3-46*pf^2+2.4e+002*pf+5.7e+002;%enthalpy of saturated feed kj/kg
hs=zeros(1,n);
hs(1,1)=0.059*p(1,1)^3-2.3*p(1,1)^2+10*p(1,1)+2.8e+003; %specific enthalpy of steam (kJ/kg)
hw=zeros(1,n);
hw(1,1)=0.24*p(1,1)^3-7.5*p(1,1)^2+1.2e+002*p(1,1)+7.1e+002;%specific enthalpy of water (kJ/kg)
Ps=zeros(1,n);
Ps(1,1)=-0.0029*p(1,1)^6+0.11*p(1,1)^5-1.8*p(1,1)^4+15*p(1,1)^3-64*p(1,1)^2+1.5e+002*p(1,1)-1.3e+002; % density of steam(kg/m^3)
Pw=zeros(1,n);
Pw(1,1)=-0.056*p(1,1)^3+1.6*p(1,1)^2-32*p(1,1)+9e+002; %density of water(kg/m^3)
U=zeros(1,n);
U(1,1)=-(6*0.0029)*p(1,1)^5+(0.11*5)*p(1,1)^4-(1.8*4)*p(1,1)^3+(15*3)*p(1,1)^2-(64*2)*p(1,1)+002; %dPs/dp
V=zeros(1,n);
V(1,1)=-(0.056*3)*p(1,1)^2+(1.6*2)*p(1,1)-32; %dPw/dp
W=zeros(1,n);
W(1,1)=(0.059*3)*p(1,1)^2-(2.3*2)*p(1,1)+10; %dhs/dp
X=zeros(1,n);
X(1,1)=(0.24*3)*p(1,1)^2-15*p(1,1)+002; %dhw/dp
ts=zeros(1,n);
ts(1,1)=0.06*p(1,1)^3-1.9*p(1,1)^2+27*p(1,1)+1.7e+002; %temperature of steam(degC)
O=zeros(1,n);
O(1,1)=(0.06*3)*p(1,1)^2-(1.9*2)*p(1,1)+27; %dts/dp
g=9.81; %acceleration due to gravity m/s^2
%first system is at equilibrium
qs=qf;
Q=qs*hs(1,1)-qf*hf;%in kW
Vsd0=5.95569;
hc=zeros(1,n);
hc(1,1)=hs(1,1)-hw(1,1);%enthalpy of condensation in kJ/kg
Vsd=zeros(1,n);
Vsd(1,1)=Vsd0-((Td*(hw(1,1)-hf)*qf)/(Ps(1,1)*hc(1,1)));%volume of steam in the drum below water level
Lr=11.739;%length of riser in m
Ldc=11.739;%length of downcomer in m
Adc=Vdc/Ldc;% area of downcomer in m^2
Ar=Vr/Lr;% area of riser in m^2
ar=zeros(1,n);
ar(1,1)=0.051;
aar=zeros(1,n);
aar(1,1)=0.1262;
neta=zeros(1,n);
neta(1,1)=ar(1,1)*(Pw(1,1)-Ps(1,1))/Ps(1,1);
daarp=zeros(1,n);
daarp(1,1)=(1/(Pw(1,1)-Ps(1,1))^2)*(Pw(1,1)*U(1,1)-Ps(1,1)*V(1,1))*(1+(Pw(1,1)/(Ps(1,1)*(1+neta(1,1)))))-(((Ps(1,1)+Pw(1,1))*log(1+neta(1,1)))/(neta(1,1)*Ps(1,1)));
daarar=zeros(1,n);
daarar(1,1)=(Pw(1,1)/(Ps(1,1)*neta(1,1)))*((log(1+neta(1,1))/neta(1,1))-(1/1+neta(1,1)));
Vwd=zeros(1,n);
Vwd(1,1)=Vwt(1,1)-Vdc-((1-aar(1,1))*Vr);
lw=zeros(1,n);
lw(1,1)=Vwd(1,1)/Ad;
ls=zeros(1,n);
ls(1,1)=Vsd(1,1)/Ad;
l=zeros(1,n);
l(1,1)=lw(1,1)+ls(1,1);
qsd=zeros(1,n);
qsd(1,1)=(Ps(1,1)*Vsd0/Td);
qdc=zeros(1,n);
qdc(1,1)=((2*Pw(1,1)*Adc*(Pw(1,1)-Ps(1,1))*g*aar(1,1)*Vr)/k)^0.5;
qct=zeros(1,n);
qct(1,1)=(hw(1,1)-hf(1,1))*qf/hc(1,1);
qr=zeros(1,n);
qr(1,1)=qdc(1,1);
e11=zeros(1,n);
e11(1,1)=Pw(1,1)-Ps(1,1);
e12=zeros(1,n);
e12(1,1)=(Vwt(1,1)*V(1,1))+(Vst(1,1)*U(1,1));
e21=zeros(1,n);
e21(1,1)=(Pw(1,1)*hw(1,1))-(Ps(1,1)*hs(1,1));
e22=zeros(1,n);
e22(1,1)=(Vwt(1,1)*(hw(1,1)*V(1,1)+Pw(1,1)*X(1,1)))+(Vst(1,1)*(hs(1,1)*U(1,1)+Ps(1,1)*W(1,1)))-Vt+(mt*Cp*O(1,1));
e32=zeros(1,n);
e32(1,1)=((Pw(1,1)*X(1,1)-ar(1,1)*hc(1,1)*V(1,1))*(1-aar(1,1))*Vr)+((((1-ar(1,1))*hc(1,1)*U(1,1))+Ps(1,1)*W(1,1))*aar(1,1)*Vr)+((Ps(1,1)+(Pw(1,1)-Ps(1,1))*ar(1,1))*hc(1,1)*Vr*daarp(1,1))-Vr+(mr*Cp*O(1,1));
e33=zeros(1,n);
e33(1,1)=((1-ar(1,1))*Ps(1,1)+ar(1,1)*Pw(1,1))*hc(1,1)*Vr*daarar(1,1);
e42=zeros(1,n);
e42(1,1)=(Vsd(1,1)*U(1,1))+((1/hc(1,1))*(Ps(1,1)*Vsd(1,1)*W(1,1)+Pw(1,1)*Vwd(1,1)*X(1,1)-Vsd(1,1)-Vwd(1,1)+md*Cp*O(1,1)))+(ar(1,1)*(1+B)*Vr*(aar(1,1)*U(1,1)+(1-aar(1,1))*V(1,1)+(Ps(1,1)-Pw(1,1))*daarp(1,1)));
e43=zeros(1,n);
e43(1,1)=ar(1,1)*(1+B)*(Ps(1,1)-Pw(1,1))*Vr*daarar(1,1);
e44=zeros(1,n);
e44(1,1)=Ps(1,1);
%give step input
qs=qs+10;
%Q=Q+10000;
%solving the dynamic equations
for i=1:n;
    dt=10;
    param=[e11(1,i),Vwt(1,i),e12(1,i),p(1,i),qf,qs,dt,e21(1,i),e22(1,i),Q,hf,hs(1,i),e32(1,i),e33(1,i),ar(1,i),...
        hc(1,i),qdc(1,i),e42(1,i),e43(1,i),e44(1,i),Vsd(1,i),Ps(1,i),Vsd0,Td,hw(1,i)];
    x0=[Vwt(1,i); p(1,i);ar(1,i); Vsd(1,i)];
    f=@(x) myfun(x,param);
    x=fsolve(f,x0);
    x(1,2)=p(1,i+1);
    x(1,3)=ar(1,i+1);
    x(1,4)=Vsd(1,i+1);
    t(1,i+1)=t(1,i)+dt;
    Vst(1,i+1)=Vt-Vwt(1,i+1);
    hs(1,i+1)=0.059*p(1,i+1)^3-2.3*p(1,i+1)^2+10*p(1,i+1)+2.8e+003; %specific enthalpy of steam (kJ/kg)
    hw(1,i+1)=0.24*p(1,i+1)^3-7.5*p(1,i+1)^2+1.2e+002*p(1,i+1)+7.1e+002;%specific enthalpy of water (kJ/kg)
    Ps(1,i+1)=-0.0029*p(1,i+1)^6+0.11*p(1,i+1)^5-1.8*p(1,i+1)^4+15*p(1,i+1)^3-64*p(1,i+1)^2+1.5e+002*p(1,i+1)-1.3e+002; % density of steam(kg/m^3)
    Pw(1,i+1)=-0.056*p(1,i+1)^3+1.6*p(1,i+1)^2-32*p(1,i+1)+9e+002; %density of water(kg/m^3)
    U(1,i+1)=-(6*0.0029)*p(1,i+1)^5+(0.11*5)*p(1,i+1)^4-(1.8*4)*p(1,i+1)^3+(15*3)*p(1,i+1)^2-(64*2)*p(1,i+1)+002; %dPs/dp
    V(1,i+1)=-(0.056*3)*p(1,i+1)^2+(1.6*2)*p(1,i+1)-32; %dPw/dp
    W(1,i+1)=(0.059*3)*p(1,i+1)^2-(2.3*2)*p(1,i+1)+10; %dhs/dp
    X(1,i+1)=(0.24*3)*p(1,i+1)^2-15*p(1,i+1)+002; %dhw/dp
    ts(1,i+1)=0.06*p(1,i+1)^3-1.9*p(1,i+1)^2+27*p(1,i+1)+1.7e+002; %temperature of steam(degC)
    O(1,i+1)=(0.06*3)*p(1,i+1)^2-(1.9*2)*p(1,i+1)+27;
    hc(1,i+1)=hs(1,i+1)-hw(1,i+1);%enthalpy of condensation in kJ/kg
    Vsd(1,i+1)=Vsd0-((Td*(hw(1,i+1)-hf)*qf)/(Ps(1,i+1)*hc(1,i+1)));
    aar(1,i+1)=(Pw(1,i+1)/(Pw(1,i+1)-Ps(1,i+1)))*(1-(Ps(1,i+1)/((Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1)))*log(1+((Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1))/Ps(1,i+1)));
    neta(1,i+1)=ar(1,i+1)*(Pw(1,i+1)-Ps(1,i+1))/Ps(1,i+1);
    daarp(1,i+1)=(1/(Pw(1,i+1)-Ps(1,i+1))^2)*(Pw(1,i+1)*U(1,i+1)-Ps(1,i+1)*V(1,i+1))*(1+(Pw(1,i+1)/(Ps(1,i+1)*(1+neta(1,i+1)))))-(((Ps(1,i+1)+Pw(1,i+1))*log(1+neta(1,i+1)))/(neta(1,i+1)*Ps(1,i+1)));
    daarar(1,i+1)=(Pw(1,i+1)/(Ps(1,i+1)*neta(1,i+1)))*((log(1+neta(1,i+1))/neta(1,i+1))-(1/1+neta(1,i+1)));
    Vwd(1,i+1)=Vwt(1,i+1)-Vdc-((1-aar(1,i+1))*Vr);
    lw(1,i+1)=Vwd(1,i+1)/Ad;
    ls(1,i+1)=Vsd(1,i+1)/Ad;
    l(1,i+1)=lw(1,i+1)+ls(1,i+1);
    qsd(1,i+1)=(Ps(1,i+1)*Vsd0/Td);
    qdc(1,i+1)=((2*Pw(1,i+1)*Adc*(Pw(1,i+1)-Ps(1,i+1))*g*aar(1,i+1)*Vr)/k)^0.5;
    qct(1,i+1)=(hw(1,i+1)-hf(1,i+1))*qf/hc(1,i+1)+(1/hc(1,i+1)*dt)*(Ps(1,i+1)*Vst(1,i+1)*W+Pw(1,i+1)*Vwt(1,i+1)*X(1,i+1)-Vt+mt*Cp*O(1,i+1))*(p(1,i+1)-p(1,i));
    qr(1,1)=qdc(1,1)-(Vr/dt)*(aar*U+(1-aar)*V+(Pw-Ps)*daarp)*(p(1,i+1)-p(1,i))+((Pw-Ps)*Vr*daarar*(ar(1,i+1)-ar(1,i+1))/dt);
    e11(1,i+1)=Pw(1,i+1)-Ps(1,i+1);
    e12(1,i+1)=(Vwt(1,i+1)*V(1,i+1))+(Vst(1,i+1)*U(1,i+1));
    e21(1,i+1)=(Pw(1,i+1)*hw(1,i+1))-(Ps(1,i+1)*hs(1,i+1));
    e22(1,i+1)=(Vwt(1,i+1)*(hw(1,i+1)*V(1,i+1)+Pw(1,i+1)*X(1,i+1)))+(Vst(1,i+1)*(hs(1,i+1)*U(1,i+1)+Ps(1,i+1)*W(1,i+1)))-Vt+(mt*Cp*O(1,i+1));
    e32(1,i+1)=((Pw(1,i+1)*X(1,i+1)-ar(1,i+1)*hc(1,i+1)*V(1,i+1))*(1-aar(1,i+1))*Vr)+((((1-ar(1,i+1))*hc(1,i+1)*U(1,i+1))+Ps(1,i+1)*W(1,i+1))*aar(1,i+1)*Vr)+((Ps(1,i+1)+(Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1))*hc(1,i+1)*Vr*daarp(1,i+1))-Vr+(mr*Cp*O(1,i+1));
    e33(1,i+1)=((1-ar(1,i+1))*Ps(1,i+1)+ar(1,i+1)*Pw(1,i+1))*hc(1,i+1)*Vr*daarar(1,i+1);
    e42(1,i+1)=(Vsd(1,i+1)*U(1,i+1))+((1/hc(1,i+1))*(Ps(1,i+1)*Vsd(1,i+1)*W(1,i+1)+Pw(1,i+1)*Vwd(1,i+1)*X(1,i+1)-Vsd(1,i+1)-Vwd(1,i+1)+md*Cp*O(1,i+1)))+(ar(1,i+1)*(1+B)*Vr*(aar(1,i+1)*U(1,i+1)+(1-aar(1,i+1))*V(1,i+1)+(Ps(1,i+1)-Pw(1,i+1))*daarp(1,i+1)));
    e43(1,i+1)=ar(1,i+1)*(1+B)*(Ps(1,i+1)-Pw(1,i+1))*Vr*daarar(1,i+1);
    e44(1,i+1)=Ps(1,i+1);
end

my function for fsolve is myfun.m

function F=myfun(x,param)
e11(1,i)=param(1);
Vwt(1,i)=param(2);
e12(1,i)=param(3);
p(1,i)=param(4);
qf=param(5);
qs=param(6);
dt=param(7);
e21(1,i)=param(8);
e22(1,i)=param(9);
Q=param(10);
hf=param(11);
hs(1,i)=param(12);
e32(1,i)=param(13);
e33(1,i)=param(14);
ar(1,i)=param(15);
hc(1,i)=param(16);
qdc(1,i)=param(17);
e42(1,i)=param(18);
e43(1,i)=param(19);
e44(1,i)=param(20);
Vsd(1,i)=param(21);
Ps(1,i)=param(22);
Vsd0=param(23);
Td=param(24);
hw(1,i)=param(25);
F=[e11(1,i)*(x(1)-Vwt(1,i))+e12(1,i)*(x(2)-p(1,i))-(qf-qs)*dt;...
    e21(1,i)*(x(1)-Vwt(1,i))+e22(1,i)*(x(2)-p(1,i))-(Q+qf*hf-qs*hs(1,i))*dt;...
    e32(1,i)*(x(2)-p(1,i))+e33(1,i)*(x(3)-ar(1,i))-(Q-ar(1,i)*hc(1,i)*qdc(1,i))*dt;...
    e42(1,i)*(x(2)-p(1,i))+e43(1,i)*(x(3)-ar(1,i))+e44(1,i)*(x(4)-Vsd(1,i))-(Ps(1,i)*(Vsd0-Vsd(1,i))/Td)*dt-((hf-hw(1,i))/hc(1,i))*dt];
end
the error shown is
??? Subscript indices must either be real positive integers or logicals.

Error in ==> myfun at 2 e11(1,i)=param(1);

Error in ==> @(x)myfun(x,param)

Error in ==> fsolve at 254 fuser = feval(funfcn{3},x,varargin{:});

Error in ==> drumboiler at 120 x=fsolve(f,x0);

Caused by: Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Please suggest a solution

Answer 1

[Mischa Kim]

BITS, in

 function F = myfun(x,param)
 e11(1,i) = param(1);

you don't have a value assigned to i. Briefly scanning through that function it looks like you can get rid of all the (1,i), for example, replace e11(1,i) by e11.

[BITS]

If i get rid of all(1,i)...will it take the e11(1,1) value for 1st time and then e11(1,2) value for second time and so on...??? and i is assigned in for =1:100

[Mischa Kim]

All the (1,i) in the function myfun, not in the main program.

[BITS]

Thank you..it helped

Answer 2

[Image Analyst]

For this line

Vwt(1,i+1)=x(1);

what is the value of (i+1) there? Is it zero or negative? Or a floating point/fractional number? It can't be.

[BITS]

since for the first iteration it will be Vwt(1,2) it will be zero

[BITS]

I have edited the question ..pls check once

[Image Analyst]

Sounds like Mischa helped you with the edited version and you got it all working. So mark his answer "Accepted" to close it out.