How to solve interdependent differtional equations

Hi

I have a system of three differetional equations

I defined my function as follow

function dydx = StratNoObs(x, y)    
D=1; %depth [m]
rho_bs=2000;%[kg/m3]
rho_ts=1000;%[kg/m3]
drho_dy=rho_bs-rho_ts/D;
rho_t = rho_ts/(rho_bs*drho_dy); % 
rho_b = rho_bs/(rho_bs*drho_dy); % [kg/m3]
U=1;
Ri = (9.81*D^2*drho_dy)/(U^2*rho_bs); 
Re = 100; 
diffusifty = 0.000025; %[m2/s]
mu =  0.001002; % [kg/(m.s)]
Sc= (mu/(rho_bs*diffusifty)); %the ratio of the coefficient of diffusion to the kinematic viscosity
d_rho=rho_t/rho_b;
b = y(1);
f = y(2);
e = y(3);

dbdx = (280.*e+(1400/3).*f)./(Ri*Re);

fA = -(2.*f)./(Re);
fB = -Ri*(-4/735).*dbdx;
fC = (1/Re).*(-(1/3).*e-(2/9).*f);
fD = -((13/11340).*e+(31/17010).*f+(1/90));
fF = -dbdx*((-21*.e-41*.f-1404)./(882.*(rho_b-rho_t+(1/42).*b)));
fG = -(1080.*b)./(Re.*Sc.*(42.*rho_b-42.*rho_t+(1/42).*b));
fI = -(1470.*rho_b-1470.*rho_t+41.*b)/(882.*rho_b-882.*rho_t+21.*b);
fJ = ((31/17010).*e+(76/25515).*f+(1/54));
dfdx = (fA+fB+fC+fD.*(fF+fG)./(-fD.*fI+fJ);

dedx = fI.*dfdx+fF+fG;

dydx = [dbdx; dfdx; dedx];
end

I am not getting the expected answer so is I am questioning if it is ok to call the dbdx and dfdx in the defintion dedx?

Thank you

Answer 1

[Sulaymon Eshkabilov]

Here is the corrected code. There were a few syntx errs and missing paranthesis.

x= linspace(0, pi);
y = [0 1 -1];
SOL=StratNoObs(x, y)   

SOL = 3×1

1.0e+05 * 0.0000 -4.9307 8.2301

function dydx = StratNoObs(x, y)    
D=1; %depth [m]
rho_bs=2000;%[kg/m3]
rho_ts=1000;%[kg/m3]
drho_dy=rho_bs-rho_ts/D;
rho_t = rho_ts/(rho_bs*drho_dy); % 
rho_b = rho_bs/(rho_bs*drho_dy); % [kg/m3]
U=1;
Ri = (9.81*D^2*drho_dy)/(U^2*rho_bs); 
Re = 100; 
diffusifty = 0.000025; %[m2/s]
mu =  0.001002; % [kg/(m.s)]
Sc= (mu/(rho_bs*diffusifty)); %the ratio of the coefficient of diffusion to the kinematic viscosity
d_rho=rho_t/rho_b;
b = y(1);
f = y(2);
e = y(3);
dbdx = (280.*e+(1400/3).*f)./(Ri*Re);
fA = -(2.*f)./(Re);
fB = -Ri*(-4/735).*dbdx;
fC = (1/Re).*(-(1/3).*e-(2/9).*f);
fD = -((13/11340).*e+(31/17010).*f+(1/90));
fF = -dbdx*((-21.*e-41.*f-1404)./(882.*(rho_b-rho_t+(1/42).*b)));
fG = -(1080.*b)./(Re.*Sc.*(42.*rho_b-42.*rho_t+(1/42).*b));
fI = -(1470.*rho_b-1470.*rho_t+41.*b)/(882.*rho_b-882.*rho_t+21.*b);
fJ = ((31/17010).*e+(76/25515).*f+(1/54));
dfdx = (fA+fB+fC+fD.*(fF+fG))./(-fD.*fI+fJ);
dedx = fI.*dfdx+fF+fG;
dydx = [dbdx; dfdx; dedx];
end

[Ahmed Aburakhia]

thanks for correcting my code. I am more concern with the methodolagy as I am not famialir with a system of ODEs. I took the matlab course but their ODEs are not interdependant. In my case I am using dbdx in solving for dfdx and using both dbdx and dfdx in solving for dedx. Is this the correct way?

[Sulaymon Eshkabilov]

Just looking at only math operations and manipulated variables here, all steps are legit. But it is said without seeing your original diff equations.