I replaced this line
with:
ind1=zeros;
ind2=zeros;
if(~isempty(ph1)&&~isempty(ph2))
for t=1:length(ph1)
ind1(t)=sub2ind(size(y),row_v(j),col_v(j),i,ph1(t));
end
for t=1:length(ph2)
ind2(t)=sub2ind(size(y),row_v(j),col_v(j),ph2(t),i);
end
sum(y(ind1(:)))-sum(y(ind2(:)))== v_link_capacity(row_v(j),col_v(j))*(x(row_v(j),i)-x(col_v(j),i));
end
if (isempty(ph1)&&~isempty(ph2))
for t=1:length(ph2)
ind2(t)=sub2ind(size(y),row_v(j),col_v(j),ph2(t),i);
end
-sum(y(ind2(:)))- v_link_capacity(row_v(j),col_v(j))*(x(row_v(j),i)-x(col_v(j),i))==0;
end
if (isempty(ph2)&&~isempty(ph1))
for t=1:length(ph1)
ind1(t)=sub2ind(size(y),row_v(j),col_v(j),i,ph1(t));
end
sum(y(ind1(:)))- v_link_capacity(row_v(j),col_v(j))*(x(row_v(j),i)-x(col_v(j),i))==0;
end