MATLAB: HOW TO SOLVE A SYSTEM OF FIRST ORDER ODE WITH ODE45 SOLVER (knowing that it contains if statment and other parameter depending on the solutions of the system)

if ode45if statment odesode's

HELLO
I would I would be grateful if you could help me.
i have a system of ODE that i want to solve with ode45;
i used if statment but it didn't seem to work properly (it gives the solution for just one condition)
the function is : function F = ODEfun(t,y)
the system is given below
i defined the ode45 solver: 
[t, y] = ode45(@ODEfun,tspan,y0);
F = [dydt(1);dydt(2);dydt(3)];
and this is the system
%===========================================================================================
rho_w = a_1 * y(2)^2 + b_1 * y(2) + c_1;
%=========================================================================================
hC = a_2*y(1) + b_2 + c_2 + d_2;
% ================================  Differential equations =========================================
% ================================  water content =========================================
if  y(2)<T_b
    dydt(1) = Q*(phi*V_total-y(1))/(phi_eff*V_total); 
elseif y(2)>= T_b
    dydt(1) = -(h*A_s/(rho_w * L_v))*(y(2) - T_b);      
elseif y(2)<T_b 
    dydt(1) = 0;                                    
end
%========================================== temperature===================
if  y(2) < T_b
    dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl  / d_inf ))...
        - Q * ( rho_w * c_w * ( y(2) - T_inf)/hC)* ((phi * V_total - y(1) ) /...
        phi_eff*V_total);
elseif  y(2)>= T_b 
    dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl  / d_inf ))...
        - ( rho_w * L_v * h * A_cyl * (y(2) - T_b) / (hC * rho_w * L_v));
elseif  y(2)<T_b  
    dydt(2) = ( 1 / hC ) * ( Pw - ( k * (y(2) - T_inf ) * A_cyl  / d_inf ));
end
%=============================================Distance ======================
if y(2)<T_b
    dydt(3) = 0;
elseif y(2)>= T_b   
    dydt(3) = (h * A_s / (rho_w * L_v * A ) ) * ( y(2) - T_b );
elseif   y(2)<T_b    
    dydt(3) = -Q / A;
end

Best Answer

Alan Stevens pointed out the conflict in your conditions. However, the only time that y(2) < T_b and y(2)>= T_b can both be false is if y(2) is nan, so your third condition would not have been reached anyhow.
However:
The mathematics of Runge Kutta assumes that the derivatives you give have continuous first and second derivatives. When you use if inside your functions, it is seldom the case that the first and second derivatives of each returned value is continuous. You would have had to have arranged your
if  y(2)<T_b
    dydt(1) = Q*(phi*V_total-y(1))/(phi_eff*V_total); 
elseif y(2)>= T_b
    dydt(1) = -(h*A_s/(rho_w * L_v))*(y(2) - T_b);      
elseif y(2)<T_b 
    dydt(1) = 0;                                    
end
to have matched the first and second derivatives of Q*(phi*V_total-y(1))/(phi_eff*V_total) and -(h*A_s/(rho_w * L_v))*(y(2) - T_b) at T == T_b, and match 0 at whatever the third condition is.
It isn't impossible... it is just messy. It can happen in some cases involving cubic spline interpolation though.
If you cannot match two derivatives your expression branches at the condition breakpoints, then ode45() will not be able to calculate the ODE properly. It will either notice the problem and complain about a singularlity, or else it will not notice the problem and will silently return the wrong value, which is a much worse circumstance (because it misleads you into thinking your answer is not nonsense when it is nonsense.)
If you cannot match the derivatives then you need to use event functions that signal termination when the conditions are crossed, and then you have to restart integration from where you left off.. this time on the other side of the boundary.
Related Question