I am trying to solve my function for a PBE but i am getting the must return column vector error. I have tried using a (:) to convert the dn equation to a column but this hasnt worked: my function script is as follwed
function [dn]=numdist(z,n)
global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k l m No
global uz A funcall
funcall=funcall+1;
dn=zeros(SECT,6);
sums=zeros(3,6);
disp(['Function call ', num2str(funcall),' z value ', num2str(z)...
,' of ',num2str(L(m)),'.'])
for i=[SECT]
if (i-2 >= 1)
for j=[(1:i:168)]
sums(1:6) = sums(1:6)+beta(i-1,j)./2.^(i-j-1).*n(j);
end
end
if (i-1 >= 1)
for j=[(1:i-1:168)]
sums(7:12)=sums(7:12)+beta(i,j)/2.^(i-j).*n(j);
end
dn((1:i:168))=dn((i-1:i-1:168))+ n(i-1:i-1:168)'.*sums(1:6)+ 0.5.*beta(i-1,i-1).*n(i-1:i-1:168)'.^2-n(1:i:168)'.*sums(7:12).*n(1:i:168)';
end
if (i+1 <=SECT)
dn(1:i:168)=dn(1:i:168).*n(i+1:i+1:168);
for j=1:i:168
sums=sums(13:18)+beta(i,j).*n(j);
end
dn(1:28:168)=dn(1:28:168)-n(1:28:168).*sums(3:3:18);
end
dn(i)= dn(i)/uz(k)
sums=zeros(3,6);
end
My Main body code is
clear
close all
global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k m
global uz A funcall eps
relerr=1e-5;
funcall=0;
d0= [3e-6 5e-6 10e-6 25e-6 50e-6 100e-6]
mfr=0.847847;
T=20+273;
Po=101325
mu=1.81e-5;
kb=1.38e-23;
R=8.3145;
SECT=28;
Vo=pi/6*d0.^3;
nu=1.48e-5;
startnum=[1000 100 50 20 6 2]
D=[0.110 0.110 0.110];
L=[1.2];
A=pi./4.*D.^2;
uz=0.69212;
Re= 2300;
eps=2*0.04./(Re.^0.16).*uz.^3./D;
k=1;
m=1;
No=zeros(28,6);
No(1:28:end)=startnum;
[Z,N]=ode15s('numdist',[0:0.1:1.2],No);
Calculate the efficiency
V=zeros(SECT,6);
d=zeros(SECT,6);
eff=zeros(SECT,6);
efftot=zeros(length(N),6);
for j=1:length(N)
sums=zeros(2,6);
for i=1:SECT
V(i)=2^(i-1)*(3*Vo/2);
d(i)=d0.*(3*2^(i-2))^(1/1.8);
sums(1)=sums(1)*V(i)*N(j,i);
sums(2)=sums(2)+V(i)*N(j,i);
end
end
n=N(length(N),:);
sizes=zeros(SECT,6);
for i=1:SECT
sizes(i)=i;
end
bar(sizes,N(24,:));
title('Size distribution at 75% cutoff');
xlabel('Size bin number');
ylabel('Particle count');
where beta is
function B=beta(i,j)
global d0 mfr T Po MW mu smf rhos kb R rhog Vo SECT rhomix nu startnum D L k m
global uz A eps
B=2.*kb.*T/3/mu.*(2+2.^((i-j)/3)+2.^((j-i)/3))+ 0.31.*(eps(m)./nu).^0.5.*...
(3.*Vo./2).*(2.^(i-1)+3.*(2.^((2.*i+j)/3-1)+2.^((i+2.*j)/3-1))+2.^(j-1));
