n=50; a=0.011;T=100; neu=0.005;Re=1e10;p0=1.013e5;r0=0;rho=1.225;
r=linspace(0.01,5,n+1);z=linspace(0,2,n+1);theta=linspace(0,2*pi,n+1);
dr=(r(n+1)-r(1))/(n+1);dz=1/(n+1);dt=0.01/(n+1);
u=zeros(n+1,n+1); v=zeros(n+1,n+1); w=zeros(n+1,n+1);p=zeros(n+1,n+1);
p(:,:)=p(:,:)+p0;
fun=@(x)(((T./(2*pi*x)).*(1-exp(-a.*(x.^2)./(2*neu)))).^2)./x;
for i=1:n+1, for j=1:n+1 % velocity defined
u(i,j)=-a.*r(i)+6.*neu.*(1-exp(-a.*(r(i).^2)./(2*neu)))./r(i);
v(i,j)=(((T./(2*pi*r(i))).*(1-exp(-a.*(r(i).^2)./(2*neu)))).^2);
w(i,j)=2*a.*z(j).*(1-3.*exp(-a.*(r(i).^2)./(2*neu)));
end, end
for i=1:n+1, for j=1:n+1 % pressure defined
p(i,j) = p(i,j)+ integral(fun,r0,r(i))-((rho*a*a/2).*(r(i).*r(i)+4.*z(j).*z(j)))-(18*rho*neu*neu./(r(i).^2)).*((1-exp(-a.*(r(i).^2)./(2*neu))).^2);
end, end
for i=1:n+1, for j=1:n+1, for k=1:n+1
U(i,j,k)=u(i,k);
V(i,j,k)=v(i,k);
W(i,j,k)=w(i,k);
P(i,j,k)=p(i,k);
end,end,end