How can i remove singular Jacobin error from this code?
william()
function william
clc
clear all format long
w=1; %wlmsn
Ha=0.1;
Re=0.4;
angle_degrees=0.1; %omega N,a,Ha,Re,angle_degrees,E,Pr
E=0.1; %Ec
Pr=0.2; %0.01
phi1=0.001;%0.02,0.03
phi2=0.001;%0.03,0.2
rof=997.1;
ro1=8933;
ro2=3970;
rocpf=4179;
rocp1= 385;
rocp2= 765;
kf=0.613;
k1=401;
k2=40
sigmaf=0.05;
sigma1=5.96*10^7;
sigma2=3.69*10^7;
A1=(phi2*rocp2)/rocpf+ (1-phi2)*((1-phi1)+(rocp1*phi1)/rocpf);
A2=((1-phi2)*((1-phi1)+(ro1*phi1)/rof)+(ro2*phi2)/rof);
A= ((sigma1*(1-2*(phi1)) + 2*sigmaf*(1-phi1))/(sigma1*(1-phi1) + 2*sigmaf*(1+phi1))); %sigmaf
A3=((sigma2*(1+2*(phi2))+2*(1-phi2)*A)/(sigma2*(1-phi2)+(2+phi2)*A))*A;
%A3=((sigma2+2*(1-phi2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1))+2*phi2*sigma2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1)))/(sigma2+2*(1+phi2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1))-phi2*sigma2);
A4=((1-phi1)^-2.5)*((1-phi2)^-2.5); %meuf
A5=(((2 * (1 + phi1) * (2 * kf – 2 * phi2 * (k2 – kf) + k2)) / (2 * kf + phi2 * (k2 – kf) + k2) + ((1 – 2 * phi1) * k1) / kf) / (2 – phi1)) * ((2 * kf – 2 * phi2 * (k2 – kf) + k2) / (2 * kf + phi2 * (k2 – kf) + k2)) + ((1 + phi1) * k1) / kf;
infinity=1;
solinit = bvpinit(linspace(0,infinity,20),zeros(5,1));
options=bvpset(‘Stats’,’on’,’RelTol’,1e-6);
lines={‘k’,’r’,’g’,’m’,’k–‘,’r–‘,’g–‘,’b–‘,’y’,’m’,’k’,’b’};
Ha=0.1;
sol1 = bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.3;
sol2 =bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.5;
sol3 =bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.7;
sol4= bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
figure(1);
plot(sol1.x,sol1.y(1,:),lines{1},sol2.x,sol2.y(1,:),lines{2},sol3.x,sol3.y(1,:),lines{3},sol4.x,sol4.y(1,:),lines{4},’LineWidth’,3,’MarkerSize’,16,’linestyle’,’-‘);
xlabel(‘zeta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(2);
plot(sol1.x,sol1.y(4,:),lines{1},sol2.x,sol2.y(4,:),lines{2},sol3.x,sol3.y(4,:),lines{3},sol4.x,sol4.y(4,:),lines{4},’LineWidth’,3,’MarkerSize’,16,’linestyle’,’-‘);
xlabel(‘zeta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘kappa’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
y1 = sol1.y;
y2 = sol2.y;
function dydx = williamode(x,y,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2)
dydx =zeros(5,1);
dydx = [y(2)
y(3)
%-(sqrt((2*y(1)*angle_degrees)^2 + y(2)^2) / (sqrt((2*y(1)*angle_degrees)^2 + y(2)^2) + w * (2 * y(2)^2 + 4 * y(1)^2 * angle_degrees^2))) * (2 * angle_degrees * Re * A2 * (1 – phi1)^2.5 * (1 – phi2)^2.5 * y(1) * y(2)+(4-A3*(1-phi1)^2.5*(1-phi2)^2.5*(Ha)^2)*y(2)*angle_degrees^2+(w*(sqrt(2*y(1)*angle_degrees)^2+y(2)^2)*((16*y(1)^3*angle_degrees^4*+8*y(1)*y(2)^2*angle_degrees^2+2*y(2)^2*y(3)+4*y(1)^2*y(3)*angle_degrees^2)/(4*y(1)*y(2)*angle_degrees^2+y(2)*y(3))+(1/((2*y(1)*angle_degrees)^2*+y(2)^2))*(80*y(1)^2*y(2)*angle_degrees^4+32*y(1)*y(2)*y(3)*angle_degrees^2+8*y(2)^3*angle_degrees^2+4*y(2)*y(3)^2)));
-(sqrt((2*y(1)*angle_degrees)^2+y(2)^2)/((sqrt((2*y(1)*angle_degrees)^2+y(2)^2))+w*(2*y(2)^2+4*y(1)^2*angle_degrees^2)))*(2*angle_degrees*Re*A2*(1-phi1)^2.5*(1-phi2)^2.5*y(1)*y(2)+(4-A3*(1-phi1)^2.5*(1-phi2)^2.5*(Ha)^2)*y(2)*angle_degrees^2+((w*(sqrt(2*y(1)*angle_degrees)^2+y(2)^2))*((16*y(1)^3*angle_degrees^4*+8*y(1)*y(2)^2*angle_degrees^2+2*y(2)^2*y(3)+4*y(1)^2*y(3)*angle_degrees^2)/(4*y(1)*y(2)*angle_degrees^2+y(2)*y(3)))+((1/((2*y(1)*angle_degrees)^2*+y(2)^2))*(80*y(1)^2*y(2)*angle_degrees^4+32*y(1)*y(2)*y(3)*angle_degrees^2+8*y(2)^3*angle_degrees^2+4*y(2)*y(3)^2))))
%-(sqrt((2 * y(1) * angle_degrees)^2 + y(2)^2) /(sqrt((2 * y(1) * angle_degrees)^2 + y(2)^2) + w * (2 * y(2)^2 + 4 * y(1)^2 * angle_degrees^2)))*(2*angle_degrees *Re*A2*(1 – phi1)^2.5*(1 – phi2)^2.5 * y(1) * y(2) +(4 – A3 * (1 – phi1)^2.5 * (1 – phi2)^2.5 * (Ha)^2) * y(2) * angle_degrees^2 +((w * (sqrt(2 * y(1) * angle_degrees)^2 + y(2)^2)) *((16 * y(1)^3 * angle_degrees^4 + 8 * y(1) * y(2)^2 * angle_degrees^2 + 2 * y(2)^2 * y(3) +4 * y(1)^2 * y(3) * angle_degrees^2) /(4 * y(1) * y(2) * angle_degrees^2 + y(2) * y(3)) +((1 / ((2 * y(1) * angle_degrees)^2 + y(2)^2)) *(80 * y(1)^2 * y(2) * angle_degrees^4 + 32 * y(1) * y(2) * y(3) * angle_degrees^2+8 * y(2)^3 * angle_degrees^2 + 4 * y(2) * y(3)^2)));
y(5)
(-1/A5)*((E*Pr*A1)/(1-phi1)^2.5*(1-phi2)^2.5)*((4*y(1)^2*angle_degrees^2+y(2)^2)+w*((4*y(1)^2*angle_degrees^2+y(2)^2)^3/2))]
end
function res = williambc(ya,yb,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,~,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2)
res =[ya(1)-1;
ya(2);
ya(5);
yb(1);
yb(4)-1;];
hold on
end
grid on
endwilliam()
function william
clc
clear all format long
w=1; %wlmsn
Ha=0.1;
Re=0.4;
angle_degrees=0.1; %omega N,a,Ha,Re,angle_degrees,E,Pr
E=0.1; %Ec
Pr=0.2; %0.01
phi1=0.001;%0.02,0.03
phi2=0.001;%0.03,0.2
rof=997.1;
ro1=8933;
ro2=3970;
rocpf=4179;
rocp1= 385;
rocp2= 765;
kf=0.613;
k1=401;
k2=40
sigmaf=0.05;
sigma1=5.96*10^7;
sigma2=3.69*10^7;
A1=(phi2*rocp2)/rocpf+ (1-phi2)*((1-phi1)+(rocp1*phi1)/rocpf);
A2=((1-phi2)*((1-phi1)+(ro1*phi1)/rof)+(ro2*phi2)/rof);
A= ((sigma1*(1-2*(phi1)) + 2*sigmaf*(1-phi1))/(sigma1*(1-phi1) + 2*sigmaf*(1+phi1))); %sigmaf
A3=((sigma2*(1+2*(phi2))+2*(1-phi2)*A)/(sigma2*(1-phi2)+(2+phi2)*A))*A;
%A3=((sigma2+2*(1-phi2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1))+2*phi2*sigma2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1)))/(sigma2+2*(1+phi2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1))-phi2*sigma2);
A4=((1-phi1)^-2.5)*((1-phi2)^-2.5); %meuf
A5=(((2 * (1 + phi1) * (2 * kf – 2 * phi2 * (k2 – kf) + k2)) / (2 * kf + phi2 * (k2 – kf) + k2) + ((1 – 2 * phi1) * k1) / kf) / (2 – phi1)) * ((2 * kf – 2 * phi2 * (k2 – kf) + k2) / (2 * kf + phi2 * (k2 – kf) + k2)) + ((1 + phi1) * k1) / kf;
infinity=1;
solinit = bvpinit(linspace(0,infinity,20),zeros(5,1));
options=bvpset(‘Stats’,’on’,’RelTol’,1e-6);
lines={‘k’,’r’,’g’,’m’,’k–‘,’r–‘,’g–‘,’b–‘,’y’,’m’,’k’,’b’};
Ha=0.1;
sol1 = bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.3;
sol2 =bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.5;
sol3 =bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.7;
sol4= bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
figure(1);
plot(sol1.x,sol1.y(1,:),lines{1},sol2.x,sol2.y(1,:),lines{2},sol3.x,sol3.y(1,:),lines{3},sol4.x,sol4.y(1,:),lines{4},’LineWidth’,3,’MarkerSize’,16,’linestyle’,’-‘);
xlabel(‘zeta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(2);
plot(sol1.x,sol1.y(4,:),lines{1},sol2.x,sol2.y(4,:),lines{2},sol3.x,sol3.y(4,:),lines{3},sol4.x,sol4.y(4,:),lines{4},’LineWidth’,3,’MarkerSize’,16,’linestyle’,’-‘);
xlabel(‘zeta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘kappa’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
y1 = sol1.y;
y2 = sol2.y;
function dydx = williamode(x,y,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2)
dydx =zeros(5,1);
dydx = [y(2)
y(3)
%-(sqrt((2*y(1)*angle_degrees)^2 + y(2)^2) / (sqrt((2*y(1)*angle_degrees)^2 + y(2)^2) + w * (2 * y(2)^2 + 4 * y(1)^2 * angle_degrees^2))) * (2 * angle_degrees * Re * A2 * (1 – phi1)^2.5 * (1 – phi2)^2.5 * y(1) * y(2)+(4-A3*(1-phi1)^2.5*(1-phi2)^2.5*(Ha)^2)*y(2)*angle_degrees^2+(w*(sqrt(2*y(1)*angle_degrees)^2+y(2)^2)*((16*y(1)^3*angle_degrees^4*+8*y(1)*y(2)^2*angle_degrees^2+2*y(2)^2*y(3)+4*y(1)^2*y(3)*angle_degrees^2)/(4*y(1)*y(2)*angle_degrees^2+y(2)*y(3))+(1/((2*y(1)*angle_degrees)^2*+y(2)^2))*(80*y(1)^2*y(2)*angle_degrees^4+32*y(1)*y(2)*y(3)*angle_degrees^2+8*y(2)^3*angle_degrees^2+4*y(2)*y(3)^2)));
-(sqrt((2*y(1)*angle_degrees)^2+y(2)^2)/((sqrt((2*y(1)*angle_degrees)^2+y(2)^2))+w*(2*y(2)^2+4*y(1)^2*angle_degrees^2)))*(2*angle_degrees*Re*A2*(1-phi1)^2.5*(1-phi2)^2.5*y(1)*y(2)+(4-A3*(1-phi1)^2.5*(1-phi2)^2.5*(Ha)^2)*y(2)*angle_degrees^2+((w*(sqrt(2*y(1)*angle_degrees)^2+y(2)^2))*((16*y(1)^3*angle_degrees^4*+8*y(1)*y(2)^2*angle_degrees^2+2*y(2)^2*y(3)+4*y(1)^2*y(3)*angle_degrees^2)/(4*y(1)*y(2)*angle_degrees^2+y(2)*y(3)))+((1/((2*y(1)*angle_degrees)^2*+y(2)^2))*(80*y(1)^2*y(2)*angle_degrees^4+32*y(1)*y(2)*y(3)*angle_degrees^2+8*y(2)^3*angle_degrees^2+4*y(2)*y(3)^2))))
%-(sqrt((2 * y(1) * angle_degrees)^2 + y(2)^2) /(sqrt((2 * y(1) * angle_degrees)^2 + y(2)^2) + w * (2 * y(2)^2 + 4 * y(1)^2 * angle_degrees^2)))*(2*angle_degrees *Re*A2*(1 – phi1)^2.5*(1 – phi2)^2.5 * y(1) * y(2) +(4 – A3 * (1 – phi1)^2.5 * (1 – phi2)^2.5 * (Ha)^2) * y(2) * angle_degrees^2 +((w * (sqrt(2 * y(1) * angle_degrees)^2 + y(2)^2)) *((16 * y(1)^3 * angle_degrees^4 + 8 * y(1) * y(2)^2 * angle_degrees^2 + 2 * y(2)^2 * y(3) +4 * y(1)^2 * y(3) * angle_degrees^2) /(4 * y(1) * y(2) * angle_degrees^2 + y(2) * y(3)) +((1 / ((2 * y(1) * angle_degrees)^2 + y(2)^2)) *(80 * y(1)^2 * y(2) * angle_degrees^4 + 32 * y(1) * y(2) * y(3) * angle_degrees^2+8 * y(2)^3 * angle_degrees^2 + 4 * y(2) * y(3)^2)));
y(5)
(-1/A5)*((E*Pr*A1)/(1-phi1)^2.5*(1-phi2)^2.5)*((4*y(1)^2*angle_degrees^2+y(2)^2)+w*((4*y(1)^2*angle_degrees^2+y(2)^2)^3/2))]
end
function res = williambc(ya,yb,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,~,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2)
res =[ya(1)-1;
ya(2);
ya(5);
yb(1);
yb(4)-1;];
hold on
end
grid on
end william()
function william
clc
clear all format long
w=1; %wlmsn
Ha=0.1;
Re=0.4;
angle_degrees=0.1; %omega N,a,Ha,Re,angle_degrees,E,Pr
E=0.1; %Ec
Pr=0.2; %0.01
phi1=0.001;%0.02,0.03
phi2=0.001;%0.03,0.2
rof=997.1;
ro1=8933;
ro2=3970;
rocpf=4179;
rocp1= 385;
rocp2= 765;
kf=0.613;
k1=401;
k2=40
sigmaf=0.05;
sigma1=5.96*10^7;
sigma2=3.69*10^7;
A1=(phi2*rocp2)/rocpf+ (1-phi2)*((1-phi1)+(rocp1*phi1)/rocpf);
A2=((1-phi2)*((1-phi1)+(ro1*phi1)/rof)+(ro2*phi2)/rof);
A= ((sigma1*(1-2*(phi1)) + 2*sigmaf*(1-phi1))/(sigma1*(1-phi1) + 2*sigmaf*(1+phi1))); %sigmaf
A3=((sigma2*(1+2*(phi2))+2*(1-phi2)*A)/(sigma2*(1-phi2)+(2+phi2)*A))*A;
%A3=((sigma2+2*(1-phi2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1))+2*phi2*sigma2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1)))/(sigma2+2*(1+phi2)*(sigma1 + 2*sigmaf – 2*phi1*(sigmaf – sigma1)) * sigmaf / (sigma1 + 2*sigmaf + phi1 * (sigmaf – sigma1))-phi2*sigma2);
A4=((1-phi1)^-2.5)*((1-phi2)^-2.5); %meuf
A5=(((2 * (1 + phi1) * (2 * kf – 2 * phi2 * (k2 – kf) + k2)) / (2 * kf + phi2 * (k2 – kf) + k2) + ((1 – 2 * phi1) * k1) / kf) / (2 – phi1)) * ((2 * kf – 2 * phi2 * (k2 – kf) + k2) / (2 * kf + phi2 * (k2 – kf) + k2)) + ((1 + phi1) * k1) / kf;
infinity=1;
solinit = bvpinit(linspace(0,infinity,20),zeros(5,1));
options=bvpset(‘Stats’,’on’,’RelTol’,1e-6);
lines={‘k’,’r’,’g’,’m’,’k–‘,’r–‘,’g–‘,’b–‘,’y’,’m’,’k’,’b’};
Ha=0.1;
sol1 = bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.3;
sol2 =bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.5;
sol3 =bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
Ha=0.7;
sol4= bvp4c(@williamode,@williambc,solinit,options,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2);
figure(1);
plot(sol1.x,sol1.y(1,:),lines{1},sol2.x,sol2.y(1,:),lines{2},sol3.x,sol3.y(1,:),lines{3},sol4.x,sol4.y(1,:),lines{4},’LineWidth’,3,’MarkerSize’,16,’linestyle’,’-‘);
xlabel(‘zeta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(2);
plot(sol1.x,sol1.y(4,:),lines{1},sol2.x,sol2.y(4,:),lines{2},sol3.x,sol3.y(4,:),lines{3},sol4.x,sol4.y(4,:),lines{4},’LineWidth’,3,’MarkerSize’,16,’linestyle’,’-‘);
xlabel(‘zeta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘kappa’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
y1 = sol1.y;
y2 = sol2.y;
function dydx = williamode(x,y,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,rof,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2)
dydx =zeros(5,1);
dydx = [y(2)
y(3)
%-(sqrt((2*y(1)*angle_degrees)^2 + y(2)^2) / (sqrt((2*y(1)*angle_degrees)^2 + y(2)^2) + w * (2 * y(2)^2 + 4 * y(1)^2 * angle_degrees^2))) * (2 * angle_degrees * Re * A2 * (1 – phi1)^2.5 * (1 – phi2)^2.5 * y(1) * y(2)+(4-A3*(1-phi1)^2.5*(1-phi2)^2.5*(Ha)^2)*y(2)*angle_degrees^2+(w*(sqrt(2*y(1)*angle_degrees)^2+y(2)^2)*((16*y(1)^3*angle_degrees^4*+8*y(1)*y(2)^2*angle_degrees^2+2*y(2)^2*y(3)+4*y(1)^2*y(3)*angle_degrees^2)/(4*y(1)*y(2)*angle_degrees^2+y(2)*y(3))+(1/((2*y(1)*angle_degrees)^2*+y(2)^2))*(80*y(1)^2*y(2)*angle_degrees^4+32*y(1)*y(2)*y(3)*angle_degrees^2+8*y(2)^3*angle_degrees^2+4*y(2)*y(3)^2)));
-(sqrt((2*y(1)*angle_degrees)^2+y(2)^2)/((sqrt((2*y(1)*angle_degrees)^2+y(2)^2))+w*(2*y(2)^2+4*y(1)^2*angle_degrees^2)))*(2*angle_degrees*Re*A2*(1-phi1)^2.5*(1-phi2)^2.5*y(1)*y(2)+(4-A3*(1-phi1)^2.5*(1-phi2)^2.5*(Ha)^2)*y(2)*angle_degrees^2+((w*(sqrt(2*y(1)*angle_degrees)^2+y(2)^2))*((16*y(1)^3*angle_degrees^4*+8*y(1)*y(2)^2*angle_degrees^2+2*y(2)^2*y(3)+4*y(1)^2*y(3)*angle_degrees^2)/(4*y(1)*y(2)*angle_degrees^2+y(2)*y(3)))+((1/((2*y(1)*angle_degrees)^2*+y(2)^2))*(80*y(1)^2*y(2)*angle_degrees^4+32*y(1)*y(2)*y(3)*angle_degrees^2+8*y(2)^3*angle_degrees^2+4*y(2)*y(3)^2))))
%-(sqrt((2 * y(1) * angle_degrees)^2 + y(2)^2) /(sqrt((2 * y(1) * angle_degrees)^2 + y(2)^2) + w * (2 * y(2)^2 + 4 * y(1)^2 * angle_degrees^2)))*(2*angle_degrees *Re*A2*(1 – phi1)^2.5*(1 – phi2)^2.5 * y(1) * y(2) +(4 – A3 * (1 – phi1)^2.5 * (1 – phi2)^2.5 * (Ha)^2) * y(2) * angle_degrees^2 +((w * (sqrt(2 * y(1) * angle_degrees)^2 + y(2)^2)) *((16 * y(1)^3 * angle_degrees^4 + 8 * y(1) * y(2)^2 * angle_degrees^2 + 2 * y(2)^2 * y(3) +4 * y(1)^2 * y(3) * angle_degrees^2) /(4 * y(1) * y(2) * angle_degrees^2 + y(2) * y(3)) +((1 / ((2 * y(1) * angle_degrees)^2 + y(2)^2)) *(80 * y(1)^2 * y(2) * angle_degrees^4 + 32 * y(1) * y(2) * y(3) * angle_degrees^2+8 * y(2)^3 * angle_degrees^2 + 4 * y(2) * y(3)^2)));
y(5)
(-1/A5)*((E*Pr*A1)/(1-phi1)^2.5*(1-phi2)^2.5)*((4*y(1)^2*angle_degrees^2+y(2)^2)+w*((4*y(1)^2*angle_degrees^2+y(2)^2)^3/2))]
end
function res = williambc(ya,yb,w,Ha,Re,angle_degrees,E,Pr,phi1,phi2,~,ro1,ro2,rocpf,rocp1,rocp2,kf,k1,k2,sigmaf,sigma1,sigma2)
res =[ya(1)-1;
ya(2);
ya(5);
yb(1);
yb(4)-1;];
hold on
end
grid on
end remove singular jacobin error MATLAB Answers — New Questions
