(what is Problem in Matlab algorithm) for solution of set of coupled 4th order ode and algebraic equation

Hello Everybody and Good Morning.

I have set of three equations that need to be solved simultaneously through matlab code to obtain results for six variables. It is actually system of equations with one 4th order equation and one an algebraic equation.The details about the problem description can be found as:

I have written the code for this problem by unsuccessfully linking ODE with other two equations as:

clear all

clc

%.......Define Guess values for Two algebraic equations No.2 and 3.........%

global gus_value mass_flux_save

gus_value =[0.0001, 334];

mass_flux_save=[];

%.......Initial Bounderies to ODE.........%

first_ivp =0.8*10^-9; % sigma at x=0

second_ivp = 1*10^-11; % d_sigma at at x=0

third_ivp = 1*10^4; % dd_sigma at at x=0

forth_ivp = 0.008; % ddd_sigma at at x=0

y0 =[first_ivp ;second_ivp;third_ivp;forth_ivp];

%.......Integeration limits/span of ODE.........%

xi =0;

xf =10*10^-9;

%...................ODE Solution Eq. No.1................%

options = odeset('RelTol',1e-10,'AbsTol',1e-20,'InitialStep',1e-15,'MaxStep',.1*(xf-xi),'Refine',5);

[xs,ys] = ode45(@(x,Y) forthorder_equa(x,Y),[xi xf],y0,options);

%....................Ploting Results............%

axis1 = subplot(2,1,1);

plot(axis1,xs,ys(:,1));

title(axis1,'x vs sigma');

ylabel(axis1,'sigma');

xlabel(axis1,'x');

axis2 = subplot(2,1,2);

plot(axis2,mass_flux_save(:,1),mass_flux_save(:,2),'--');

title(axis2,'x vs mass flux');

ylabel(axis2,'mass flux');

xlabel(axis2,'x');

%.......ENd.........%

function dy =forthorder_equa(x,Y)

%..............Properties Defination...............%

v = 0.000282/958.31;%viscosity

sur_ten = 7.264*10^-2;%Surface Tension

A = -8.7*10^-20;%Hammaker Constant

%..............Calling Global Varaibles...............%

global gus_value mass_flux_save

delta_sigma = Y;

%..............Calling set of Non linear Equation Function...............%

non_fun = @(u)nonlinear_equation_p_m_T(u,delta_sigma);

%..............Solving set of Non linear Equations Eq.2 and 3...............%

intial_gus =[ gus_value(1),gus_value(2)];

options = optimoptions('fsolve','Display','iter');

options = optimoptions(options,'MaxIterations', 1000);

options = optimoptions(options,'OptimalityTolerance', 1e-10);

options = optimoptions(options,'FunctionTolerance', 1e-3);

options = optimoptions(options,'StepTolerance', 1e-3);

options = optimoptions(options,'FunValCheck', 'off');

options = optimoptions(options,'Algorithm', 'trust-region-dogleg','Display','off');

i = fsolve(non_fun,intial_gus,options);

mass_flux = i(1);

gus_value =[i(1), i(2)];

mass_flux_save = [mass_flux_save;x,i(1)] ;

%..............Decompostion of 4th order ODE Equation No.1...............%

first=Y(2);

second=Y(3);

third=Y(4);

forth=-(3*(A*(Y(2)^2 + 1)^(7/2)*Y(2)^2 + 5*sur_ten*Y(1)^5*Y(2)^2*Y(3)^3 - mass_flux*v*(Y(2)^2 + 1)^(7/2)*Y(1)^2 + sur_ten*Y(1)^4*Y(2)*Y(4) - sur_ten*(Y(2)^2 + 1)*Y(1)^5*Y(3)^3 + 2*sur_ten*Y(1)^4*Y(2)^3*Y(4) + sur_ten*Y(1)^4*Y(2)^5*Y(4) - A*(Y(2)^2 + 1)^(7/2)*Y(1)*Y(3) - 3*sur_ten*(Y(2)^2 + 1)*Y(1)^4*Y(2)^2*Y(3)^2 - 3*sur_ten*(Y(2)^2 + 1)*Y(1)^5*Y(2)*Y(3)*Y(4)))/(sur_ten*Y(1)^5*(2*Y(2)^2 + Y(2)^4 + 1));

dy = [first;second;third;forth];

function F = nonlinear_equation_p_m_T(u,delta_sigma)

%..............Properties Defination...............%

M = 0.018;

R = 8.314;

P_v= 1.88 *10^4;

T_v= 333;

T_w = 334;

A = -8.7*10^-20;

k_l = .68;

h_fg = 2.358*10^6;

sur_ten = 7.264*10^-2;

Vm=0.022;

%..............Finding Constants...............%

a=((M/(2*pi*R)^0.5)*((P_v*Vm)/(R*T_v)));

b=((M/(2*pi*R)^0.5)*((P_v*M*h_fg)/(R)));

%..............Calling delta_sigma from forth order equation function...............%

P_d = A/delta_sigma(1)^3;

P_c = sur_ten*(delta_sigma(3)*(1+(delta_sigma(2))^2)^-1.5);

%.............Non Linear Equations(Eq.2 and 3)..........................%

F(1) = u(1)*(u(2)^1.5)-(a*(u(2)-T_v))+(b*(P_d+P_c));

F(2) = u(1)*h_fg-(((k_l*(T_w-u(2)))/delta_sigma(1)));

The expected Results of this program should be as following:

However, I failed to get the required results as:

Please have a kind look to my problem and help me to figure it out the mistakes in the alogrithm. If you can help me with any better code for the solution of this problem, i will be more than thankful and indebted to you.

Cheers, Together we can

[J. Alex Lee]

https://www.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html