TFV.m

x=logspace(-4,4);

for i=1:50
if x(i)<2
    Ne(i)=24/x(i);
elseif x(i) >= 2 & x(i) < 500
    Ne(i)=18.5*x(i)^(-0.6);
else    Ne(i)=0.45;
end
end

for i=1:50
    Ar(i)=Ne(i)*x(i)^2;
end

for i=1:50
    Rna(i)=Ne(i)/x(i);
end


loglog(x,Ne,'r')

hold on

loglog(x,Ar,'g')

loglog(x,Rna,'b')

axis([1e-4 1e4 1e-5 1e7]);

set(gca,'XGrid','on','YGrid','on','XMinorGrid','on','YMinorGrid','on', ...
    'XMinorTick','on','YMinorTick','on',...
    'XTick',[1e-5  1e-4 1e-3  1e-2 1e-1 1 10 1e2 1000  1e4],...
    'YTick',[1e-5  1e-4 1e-3  1e-2 1e-1 1 10 1e2 1000  1e4  ...
    1e5 1e6 1e7])

nu=1.4e-3; 
Dp = 5e-3; 
rho2 = 1260; 
rho1 = 1042; 
g = 9.81;

clear Ar Ne x

Ar = 4/3*Dp^3*rho1*abs(rho2 - rho1)*g/nu^2

options=optimset('Display','off');

Re=fsolve(@equation,500,options)

vel=feval(@tfvfunct,Re)