Error in bessel code

4 views (last 30 days)
Subha S
Subha S on 29 Aug 2024
Commented: Subha S on 6 Sep 2024
I am trying to recreate the graphs in the paper "Effects of porosity in four-layered non-linear blood rheology in constricted narrow arteries with clinical applications" by Afiqah Wajihah S. and D.S. Sankar, specifically figure 9 in page 10. I'm attaching the paper for your reference. They have specified to find the constants C2, C4, C5 and C6 (which I have obtained) through MATLAB.
syms alpha2 Da alpha We m mu hm phi beta L0 R0 delta d h n z P0 Rc Rp Rb Rd sigma % Parameters
syms C2 C4 C5 C6
f1=C2-C4+((P0*(Rc)^2)/4)+((We^2*(m-1)*(P0)^3)/(64*(mu^2)*(hm^2)))*(1-(mu*hm*(Rc^2)));
f2=C4-(C5*besseli(0,sigma*Rp))-(C6*besselk(0,sigma*Rp))-(P0*((Rp^2)+(4*Da))/4);
f3=((beta*phi*(besseli(0,sigma*Rp))-(sigma*sqrt(Da)*(besseli(1,sigma*Rp))))*C5)+(((sigma*sqrt(Da)*(besselk(1,sigma*Rp)))+(beta*phi*(besselk(0,sigma*Rp))))*C6)-((P0*phi*(((Rp)*sqrt(Da))-2*beta*Da))/2);
f4=(((sigma*sqrt(Da)*(besseli(1,sigma*Rb)))-(alpha*(besseli(0,sigma*Rb))))*C5)-(((sigma*sqrt(Da)*(besselk(1,sigma*Rb)))+(alpha*(besselk(0,sigma*Rb))))*C6);
sol = solve([f1 f2 f3 f4],[C1 C2 C3 C4 C5 C6], 'ReturnConditions',true)
[sol.C2 sol.C4 sol.C5 sol.C6]
Then find the value of the pressure gradient by taking the total fow rate as 1.
% Parameters
alpha2=1.1;
Da=0.01; %(0.01-0.3)
alpha=0.01; %(0.01-0.1)
We=0.5; %(0.5-0.9)
m=2; %(2-6) %(5-6)
mu=0.2;
hm=0.2; %(0.2-0.5)
phi=0.5;
beta=0.1;
L0=1;
R0=1;
delta=0.015;
d=4;
h=1.05;
n=2;
s=(delta/((R0)*((L0)^n)))*(n^(n/(n-1)))*(1/(n-1)); z=(d+((L0)/(n^(1/(n-1)))));
R=[0.86,0.90,0.95,1];
for i=1:length(R)
R(i)=((R0)*(1-(s*((((L0)^(n-1))*(z-d))-((z-d)^n)))));
end
syms P0
sigma=1/((alpha2*Da)^(1/2));
C2=((P0^3*We^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*We^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*We^2*m*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (P0^3*We^2*m*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (Da^(1/2)*P0^3*We^2*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*We^2*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (64*Da*P0*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (64*Da*P0*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (16*P0*(R(1))^2*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (16*P0*(R(1))^2*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (16*P0*(R(2))^2*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (16*P0*(R(2))^2*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - P0^3*We^2*alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0^3*We^2*alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (P0^3*(R(1))^2*We^2*hm*mu*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (P0^3*(R(1))^2*We^2*hm*mu*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + P0^3*We^2*alpha*beta*m*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - P0^3*We^2*alpha*beta*m*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*We^2*alpha*m*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*alpha*m*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (64*Da^(3/2)*P0*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (64*Da^(3/2)*P0*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + 16*P0*(R(1))^2*alpha*beta*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - 16*P0*(R(1))^2*alpha*beta*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - 16*P0*(R(2))^2*alpha*beta*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 16*P0*(R(2))^2*alpha*beta*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (16*Da^(1/2)*P0*(R(1))^2*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(1))^2*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(2))^2*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(2))^2*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (P0^3*(R(1))^2*We^2*hm*m*mu*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*(R(1))^2*We^2*hm*m*mu*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (32*Da*P0*(R(2))*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (32*Da*P0*(R(2))*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - 32*Da^(1/2)*P0*(R(2))*alpha*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 32*Da^(1/2)*P0*(R(2))*alpha*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*P0^3*We^2*beta*m*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*We^2*beta*m*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + P0^3*(R(1))^2*We^2*alpha*beta*hm*mu*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - P0^3*(R(1))^2*We^2*alpha*beta*hm*mu*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*mu*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*mu*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(1))^2*beta*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(1))^2*beta*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(2))^2*beta*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(2))^2*beta*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*m*mu*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*m*mu*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*mu*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*mu*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - P0^3*(R(1))^2*We^2*alpha*beta*hm*m*mu*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0^3*(R(1))^2*We^2*alpha*beta*hm*m*mu*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*m*mu*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*m*mu*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))/(64*hm^2*mu^2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
C4=((P0*(R(2))^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0*(R(2))^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (4*Da*P0*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (4*Da*P0*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (4*Da^(3/2)*P0*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (4*Da^(3/2)*P0*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (2*Da*P0*(R(2))*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (2*Da*P0*(R(2))*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - 2*Da^(1/2)*P0*(R(2))*alpha*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 2*Da^(1/2)*P0*(R(2))*alpha*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - P0*(R(2))^2*alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0*(R(2))^2*alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*P0*(R(2))^2*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0*(R(2))^2*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0*(R(2))^2*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0*(R(2))^2*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))/(4*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
C5=-(P0*phi*(alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2)) + (Da^(1/2)*besselk(1, (R(3))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))*(2*Da*beta - Da^(1/2)*(R(2))))/(2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
C6=(P0*phi*(alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2)) - (Da^(1/2)*besseli(1, (R(3))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))*(2*Da*beta - Da^(1/2)*(R(2))))/(2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
f1=(P0/(8*(mu^2)*(hm^2)))+((P0/(8*(mu^2)*(hm^2)))*(1+(mu*hm*(R(1))^2))*(log(1+(mu*hm*(R(1))^2))))-((P0/(8*(mu^2)*(hm^2)))*(1+(mu*hm*((R(1))^2))))+(((We^2)*(m-1)*((P0)^3)*log(1+(mu*hm*(R(1))^2)))/(64*(mu^3)*(hm^3)))+(((We^2)*(m-1)*((P0)^3)*(1+(mu*hm*((R(1))^2)))*mu*hm*((R(1))^2))/(128*(mu^3)*(hm^3)*(1+(mu*hm*((R(1))^2)))))+((C2*(R(1))^2)/2)-((P0*(R(2))^4)/16)+((P0*(R(1))^4)/16)+((C4*(R(2))^2)/2)-((C4*(R(1))^2)/2)+((C5*(R(3))*besseli(1,sigma*(R(3)))))-((C5*(R(2))*besseli(1,sigma*(R(2)))))-((C6*(R(3))*besselk(1,sigma*(R(3)))))+((C6*(R(2))*besselk(1,sigma*(R(2)))))-((P0*Da*((R(2))^2))/2)+((P0*Da*(((R(4)))^2))/2)-1;
sol = vpasolve(f1,P0)
[sol.P0]
Then, the velocity profile is obtained. However, the profile I'm obtaining doesn't match with the pattern given in the paper. I'm attaching the code here for your reference:
% Parameters
alpha2=1.1;
Da=0.01; %(0.01-0.3)
alpha=0.01; %(0.01-0.1)
P0=0.68357421327493338931529679176865;
We=0.5; %(0.5-0.9)
m=2; %(2-6)
mu=0.2;
hm=0.2; %(0.2-0.5)
phi=0.5;
beta=0.1;
L0=1;
R0=1;
delta=0.015;
d=4;
h=1.05;
n=2;
s=(delta/(R0)*((L0)^n))*(n^(n/(n-1)))*(1/(n-1)); z=d+((L0)/(n^(1/(n-1))));
R=[0.86,0.90,0.95,1];
for i=1:length(R)
R(i)=((R0)*(1-(s*((((L0)^(n-1))*(z-d)-((z-d)^n))))));
end
sigma=1/((alpha2*Da)^(1/2));
C2=((P0^3*We^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*We^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*We^2*m*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (P0^3*We^2*m*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (Da^(1/2)*P0^3*We^2*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*We^2*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (64*Da*P0*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (64*Da*P0*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (16*P0*(R(1))^2*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (16*P0*(R(1))^2*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (16*P0*(R(2))^2*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (16*P0*(R(2))^2*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - P0^3*We^2*alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0^3*We^2*alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (P0^3*(R(1))^2*We^2*hm*mu*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (P0^3*(R(1))^2*We^2*hm*mu*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + P0^3*We^2*alpha*beta*m*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - P0^3*We^2*alpha*beta*m*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*We^2*alpha*m*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*alpha*m*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (64*Da^(3/2)*P0*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (64*Da^(3/2)*P0*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + 16*P0*(R(1))^2*alpha*beta*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - 16*P0*(R(1))^2*alpha*beta*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - 16*P0*(R(2))^2*alpha*beta*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 16*P0*(R(2))^2*alpha*beta*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (16*Da^(1/2)*P0*(R(1))^2*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(1))^2*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(2))^2*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(2))^2*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (P0^3*(R(1))^2*We^2*hm*m*mu*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*(R(1))^2*We^2*hm*m*mu*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (32*Da*P0*(R(2))*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (32*Da*P0*(R(2))*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - 32*Da^(1/2)*P0*(R(2))*alpha*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 32*Da^(1/2)*P0*(R(2))*alpha*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*P0^3*We^2*beta*m*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*We^2*beta*m*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + P0^3*(R(1))^2*We^2*alpha*beta*hm*mu*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - P0^3*(R(1))^2*We^2*alpha*beta*hm*mu*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*mu*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*mu*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(1))^2*beta*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(1))^2*beta*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(2))^2*beta*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(2))^2*beta*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*m*mu*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*m*mu*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*mu*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*mu*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - P0^3*(R(1))^2*We^2*alpha*beta*hm*m*mu*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0^3*(R(1))^2*We^2*alpha*beta*hm*m*mu*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*m*mu*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*m*mu*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))/(64*hm^2*mu^2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
C4=((P0*(R(2))^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0*(R(2))^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (4*Da*P0*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (4*Da*P0*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (4*Da^(3/2)*P0*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (4*Da^(3/2)*P0*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (2*Da*P0*(R(2))*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (2*Da*P0*(R(2))*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - 2*Da^(1/2)*P0*(R(2))*alpha*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 2*Da^(1/2)*P0*(R(2))*alpha*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - P0*(R(2))^2*alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0*(R(2))^2*alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*P0*(R(2))^2*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0*(R(2))^2*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0*(R(2))^2*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0*(R(2))^2*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))/(4*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
C5=-(P0*phi*(alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2)) + (Da^(1/2)*besselk(1, (R(3))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))*(2*Da*beta - Da^(1/2)*(R(2))))/(2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
C6=(P0*phi*(alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2)) - (Da^(1/2)*besseli(1, (R(3))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))*(2*Da*beta - Da^(1/2)*(R(2))))/(2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
r=linspace(0,1,50);
uc=(P0*log(1+(mu*hm*(((R(1))^2)-(r.^2))))/(4*mu*hm))-((We^2)*(m-1)*(P0^3)*((1./(1+(mu*hm*(((R(1))^2)-(r.^2)))))-((1+(mu*hm*((R(1))^2)))./((2*(1+(mu*hm*(((R(1))^2)-(r.^2)))).^2))))/(32*(mu^2)*(hm^2)))+C2;
un=(-(P0*(r.^2))/4)+C4;
ub=(C5*besseli(0,sigma*r))+(C6*besselk(0,sigma*r))+(P0*Da);
ud=P0*Da;
u=uc+un+ub+ud;
plot(u,r)
xlabel('u')
ylabel('r')

Sign in to answer this question.

Accepted Answer

Divyajyoti Nayak
Divyajyoti Nayak on 2 Sep 2024
Edited: Divyajyoti Nayak on 4 Sep 2024
Hi @Subha S,
From what I understand, your plots are not matching the figure given in the research paper. This is because of a simple mistake in making your velocity vector. The velocity vector ‘u’ is a function of the radius ‘r’. So when the radius lies in the core region the velocity is uc , in the plasma region up and so on. In your code you have added all the velocities up which gives the wrong result.
Here's the code to fix the issue:
u=[uc(r<=R(1)),un(and(r<=R(2),r>R(1))),ub(and(r<=R(3),r>R(2))),ud*ones(1,sum(r>R(3)))];
And here's the result after making this change:
Another issue I noticed is in the definition of ‘sigma’. In the paper, ‘sigma’ is defined using ‘α1’ but in your code you have used ‘alpha2’. Also the definitions of the constants do not match with the result from the symbolic calculation in your first code. I changed the definition of ‘sigma’ and the constants in the code snippet below:
alpha1 = 0.01;
%Changed from alpha2 to alpha1
sigma=1/((alpha1*Da).^(1/2));
%result from Symbolic calculation
C2 = (Da*P0^3*We^2*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*P0^3*We^2*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - P0^3*We^2*alpha*beta*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + P0^3*We^2*alpha*beta*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*P0^3*We^2*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*P0^3*We^2*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) - Da*P0^3*We^2*m*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) + Da*P0^3*We^2*m*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) + 64*Da^2*P0*hm^2*mu^2*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - 64*Da^2*P0*hm^2*mu^2*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - 16*Da*P0*Rc^2*hm^2*mu^2*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) + 16*Da*P0*Rc^2*hm^2*mu^2*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) + 16*Da*P0*Rp^2*hm^2*mu^2*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - 16*Da*P0*Rp^2*hm^2*mu^2*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) + P0^3*We^2*alpha*beta*m*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) - P0^3*We^2*alpha*beta*m*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) + Da^(1/2)*P0^3*We^2*alpha*m*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) + Da^(1/2)*P0^3*We^2*alpha*m*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da^(1/2)*P0^3*We^2*beta*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*P0^3*We^2*beta*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) - 64*Da^(3/2)*P0*alpha*hm^2*mu^2*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - 64*Da^(3/2)*P0*alpha*hm^2*mu^2*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) - 32*Da^(1/2)*P0*Rp*alpha*hm^2*mu^2*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + 32*Da^(1/2)*P0*Rp*alpha*hm^2*mu^2*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*P0^3*We^2*beta*m*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) - Da^(1/2)*P0^3*We^2*beta*m*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) + 16*P0*Rc^2*alpha*beta*hm^2*mu^2*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) - 16*P0*Rc^2*alpha*beta*hm^2*mu^2*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - 16*P0*Rp^2*alpha*beta*hm^2*mu^2*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + 16*P0*Rp^2*alpha*beta*hm^2*mu^2*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da*P0^3*Rc^2*We^2*hm*mu*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) + Da*P0^3*Rc^2*We^2*hm*mu*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) + 16*Da^(1/2)*P0*Rc^2*alpha*hm^2*mu^2*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) + 16*Da^(1/2)*P0*Rc^2*alpha*hm^2*mu^2*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) - 16*Da^(1/2)*P0*Rp^2*alpha*hm^2*mu^2*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - 16*Da^(1/2)*P0*Rp^2*alpha*hm^2*mu^2*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + 32*Da*P0*Rp*hm^2*mu^2*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + 32*Da*P0*Rp*hm^2*mu^2*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) + P0^3*Rc^2*We^2*alpha*beta*hm*mu*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) - P0^3*Rc^2*We^2*alpha*beta*hm*mu*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) + Da^(1/2)*P0^3*Rc^2*We^2*alpha*hm*mu*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) + Da^(1/2)*P0^3*Rc^2*We^2*alpha*hm*mu*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da*P0^3*Rc^2*We^2*hm*m*mu*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*P0^3*Rc^2*We^2*hm*m*mu*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - 16*Da^(1/2)*P0*Rc^2*beta*hm^2*mu^2*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) - 16*Da^(1/2)*P0*Rc^2*beta*hm^2*mu^2*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) + 16*Da^(1/2)*P0*Rp^2*beta*hm^2*mu^2*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + 16*Da^(1/2)*P0*Rp^2*beta*hm^2*mu^2*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*P0^3*Rc^2*We^2*alpha*hm*m*mu*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*P0^3*Rc^2*We^2*alpha*hm*m*mu*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) - Da^(1/2)*P0^3*Rc^2*We^2*beta*hm*mu*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) - Da^(1/2)*P0^3*Rc^2*We^2*beta*hm*mu*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) - P0^3*Rc^2*We^2*alpha*beta*hm*m*mu*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + P0^3*Rc^2*We^2*alpha*beta*hm*m*mu*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) + Da^(1/2)*P0^3*Rc^2*We^2*beta*hm*m*mu*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*P0^3*Rc^2*We^2*beta*hm*m*mu*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma))/(64*hm^2*mu^2*(Da*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - alpha*beta*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + alpha*beta*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma)));
C4 = (4*Da^2*P0*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - 4*Da^2*P0*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - 4*Da^(3/2)*P0*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - 4*Da^(3/2)*P0*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da*P0*Rp^2*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*P0*Rp^2*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) + 2*Da*P0*Rp*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + 2*Da*P0*Rp*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma) - 2*Da^(1/2)*P0*Rp*alpha*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + 2*Da^(1/2)*P0*Rp*alpha*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - P0*Rp^2*alpha*beta*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + P0*Rp^2*alpha*beta*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*P0*Rp^2*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*P0*Rp^2*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da^(1/2)*P0*Rp^2*beta*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*P0*Rp^2*beta*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma))/(4*(Da*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - alpha*beta*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + alpha*beta*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma)));
C5 = -(P0*phi*(alpha*besselk(0, Rb*sigma) + Da^(1/2)*sigma*besselk(1, Rb*sigma))*(2*Da*beta - Da^(1/2)*Rp))/(2*(Da*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - alpha*beta*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + alpha*beta*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma)));
C6 = (P0*phi*(alpha*besseli(0, Rb*sigma) - Da^(1/2)*sigma*besseli(1, Rb*sigma))*(2*Da*beta - Da^(1/2)*Rp))/(2*(Da*sigma^2*besseli(1, Rb*sigma)*besselk(1, Rp*sigma) - Da*sigma^2*besselk(1, Rb*sigma)*besseli(1, Rp*sigma) - alpha*beta*phi*besseli(0, Rb*sigma)*besselk(0, Rp*sigma) + alpha*beta*phi*besselk(0, Rb*sigma)*besseli(0, Rp*sigma) - Da^(1/2)*alpha*sigma*besseli(0, Rb*sigma)*besselk(1, Rp*sigma) - Da^(1/2)*alpha*sigma*besselk(0, Rb*sigma)*besseli(1, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besseli(1, Rb*sigma)*besselk(0, Rp*sigma) + Da^(1/2)*beta*phi*sigma*besselk(1, Rb*sigma)*besseli(0, Rp*sigma)));
%Your code
% C2=((P0^3*We^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*We^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*We^2*m*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (P0^3*We^2*m*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (Da^(1/2)*P0^3*We^2*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*We^2*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (64*Da*P0*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (64*Da*P0*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (16*P0*(R(1))^2*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (16*P0*(R(1))^2*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (16*P0*(R(2))^2*hm^2*mu^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (16*P0*(R(2))^2*hm^2*mu^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - P0^3*We^2*alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0^3*We^2*alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (P0^3*(R(1))^2*We^2*hm*mu*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (P0^3*(R(1))^2*We^2*hm*mu*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + P0^3*We^2*alpha*beta*m*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - P0^3*We^2*alpha*beta*m*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*We^2*alpha*m*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*alpha*m*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*We^2*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (64*Da^(3/2)*P0*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (64*Da^(3/2)*P0*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + 16*P0*(R(1))^2*alpha*beta*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - 16*P0*(R(1))^2*alpha*beta*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - 16*P0*(R(2))^2*alpha*beta*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 16*P0*(R(2))^2*alpha*beta*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (16*Da^(1/2)*P0*(R(1))^2*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(1))^2*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(2))^2*alpha*hm^2*mu^2*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(2))^2*alpha*hm^2*mu^2*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (P0^3*(R(1))^2*We^2*hm*m*mu*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0^3*(R(1))^2*We^2*hm*m*mu*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (32*Da*P0*(R(2))*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (32*Da*P0*(R(2))*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - 32*Da^(1/2)*P0*(R(2))*alpha*hm^2*mu^2*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 32*Da^(1/2)*P0*(R(2))*alpha*hm^2*mu^2*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*P0^3*We^2*beta*m*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*We^2*beta*m*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + P0^3*(R(1))^2*We^2*alpha*beta*hm*mu*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) - P0^3*(R(1))^2*We^2*alpha*beta*hm*mu*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*mu*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*mu*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(1))^2*beta*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (16*Da^(1/2)*P0*(R(1))^2*beta*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(2))^2*beta*hm^2*mu^2*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (16*Da^(1/2)*P0*(R(2))^2*beta*hm^2*mu^2*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*m*mu*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*alpha*hm*m*mu*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*mu*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*mu*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - P0^3*(R(1))^2*We^2*alpha*beta*hm*m*mu*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0^3*(R(1))^2*We^2*alpha*beta*hm*m*mu*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) + (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*m*mu*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0^3*(R(1))^2*We^2*beta*hm*m*mu*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))/(64*hm^2*mu^2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
% C4=((P0*(R(2))^2*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (P0*(R(2))^2*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 + (4*Da*P0*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (4*Da*P0*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (4*Da^(3/2)*P0*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (4*Da^(3/2)*P0*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (2*Da*P0*(R(2))*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (2*Da*P0*(R(2))*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - 2*Da^(1/2)*P0*(R(2))*alpha*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + 2*Da^(1/2)*P0*(R(2))*alpha*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - P0*(R(2))^2*alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + P0*(R(2))^2*alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*P0*(R(2))^2*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*P0*(R(2))^2*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0*(R(2))^2*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*P0*(R(2))^2*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))/(4*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
% C5=-(P0*phi*(alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2)) + (Da^(1/2)*besselk(1, (R(3))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))*(2*Da*beta - Da^(1/2)*(R(2))))/(2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
% C6=(P0*phi*(alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2)) - (Da^(1/2)*besseli(1, (R(3))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2))*(2*Da*beta - Da^(1/2)*(R(2))))/(2*((besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - (besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/alpha2 - alpha*beta*phi*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)) + alpha*beta*phi*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)) - (Da^(1/2)*alpha*besseli(0, (R(3))/(Da*alpha2)^(1/2))*besselk(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) - (Da^(1/2)*alpha*besselk(0, (R(3))/(Da*alpha2)^(1/2))*besseli(1, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besseli(1, (R(3))/(Da*alpha2)^(1/2))*besselk(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2) + (Da^(1/2)*beta*phi*besselk(1, (R(3))/(Da*alpha2)^(1/2))*besseli(0, (R(2))/(Da*alpha2)^(1/2)))/(Da*alpha2)^(1/2)));
And here's the result:
As you can see the shape of the graph is correct but the values are wrong. This is because I took an arbitrary value of ‘α1’. You’ll have to apply the correct function for ‘α1’ as it is not given in the paper.
Hope this helps in some way.
  3 Comments
Show 1 older comment
Divyajyoti Nayak
Divyajyoti Nayak on 4 Sep 2024
Hi @Subha S,
I forgot to mention another error. Your equation for 'uc' has a slight mistake. There should be a '+' before 'We^2' according to the paper.
uc=(P0*log(1+(mu*hm*((((Rc))^2)-(r.^2))))/(4*mu*hm)) + ... % + here instead of -
((We^2)*(m-1)*(P0^3)*((1./(1+(mu*hm*((((Rc))^2)-(r.^2)))))-((1+(mu*hm*(((Rc))^2)))./((2*(1+(mu*hm*((((Rc))^2)-(r.^2)))).^2))))/(32*(mu^2)*(hm^2)))+C2;
That should give you the graphs, similar to what I obtained.
Subha S
Subha S on 6 Sep 2024
Hello Divya.
I've changed the '-' to '+' in uc equation before We^2. I'm now able to get the desired graph. Thank you so much

Sign in to comment.

More Answers (0)

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Tags

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!