У меня есть набор из трех уравнений, которые необходимо решить одновременно, чтобы получить результаты для шести переменных. На самом деле это система уравнений с одним уравнением 4-го порядка и одним алгебраическим уравнением. Подробности описания проблемы: 
Я написал код для этой проблемы, неудачно связав ODE с двумя другими уравнениями в виде:
%.......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)));
Ожидаемые результаты этой программы должны быть следующими: 
Однако я не смог получить требуемые результаты в виде: 