Must return a column vector error when using "ode45"

%% Clean the Workspace
clc
clear all

%% Format
format long

%% Variables
global m_ExtraSteam U0_shell V rho_0 U_steam

T0_shell = 105;                                       % initial shell side temp, degC
P0_shell = 1.2;                                  % initial pressure, bara
U0_shell = XSteam('uV_p', P0_shell);

Plimit  = 3.3;                                        % max pressure, bara

MW      = 0.018015;                             % kg/mol
R       = 0.00008205736;                        % bar.m3/mol.K
V       = 233.21;                                   % m3
rho_0   = XSteam('rhoV_p', P0_shell); 

tb      = 20;                                  % Upper time interval limit

m_ExtraSteam  = 35.08;           % Flowrate of extra steam due to ST trip, kg/sec

t_interval = [0.001 tb];                                     % Time interval

%% Solve the DE
[tsol,Usol] = ode45(@(t, U_steam) cfunc(t, U_steam) , t_interval , U0_shell);      % Temperature with time, degC

Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.

Equation solved at initial point.

fsolve completed because the vector of function values at the initial
point is near zero as measured by the value of the function tolerance,
and the problem appears regular as measured by the gradient.

Error using odearguments (line 94)
@(T,U_STEAM)CFUNC(T,U_STEAM) must return a column vector.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%% Dynamic Parameter Calculation
Ptube_in = 7;
Ptube_out = 6.5;

m_tube = [707.9 707.9 707.9 707.9];
T_tube_in = [91.7 91.7 91.7 91.7];

    for m=1:lenght(Usol)
        U_steam = Usol(m);
        Tsol(m) = fsolve(@TfinderwU,100);
    end

    for p=5:length(Tsol)
T_tube_in(p) = 80.1;
m_tube(p) = 1131.2;
    end

T_tube_out = 98;

        for k = 1:length(Tsol)
            deltaH_tube(k) = XSteam('h_pT', Ptube_out, T_tube_out) - XSteam('h_pT', Ptube_in, T_tube_in(k));        %  Specific enthalpy difference between the inlet and outlet of the tube side, kJ/kg
   
            if ((m_tube(k)*deltaH_tube(k))/(XSteam('hV_T', Tsol(k)) - XSteam('hL_T', Tsol(k)))) >= m_ExtraSteam
                m_ExtraSteam_Condensation(k) = m_ExtraSteam;
            else
                m_ExtraSteam_Condensation(k) = (m_tube(k)*deltaH_tube(k))/(XSteam('hV_T', Tsol(k)) - XSteam('hL_T', Tsol(k)));
            end

        end

m_steam_accumulated = (m_ExtraSteam - m_ExtraSteam_Condensation)';                 % Steam accumulated in the shell side, kg/sec

nfor(1) = m_steam_accumulated(1)*(tsol(1))/MW;

        for k1 = 2:length(Tsol)
            nfor(k1) =  m_steam_accumulated(k1)*(tsol(k1)-tsol(k1-1))/MW;
        end


%% Final Pressure Calculation



%% Plotting
figure('Name','Temperature, Pressure vs Time','NumberTitle','off');
yyaxis left % subplot(1,2,1)
plot(tsol,Tsol)
xlabel('Time (s)')
ylabel('Temperature (C)')
yyaxis right % subplot(1,2,2)
plot(tsol,Psol)
xlabel('Time (s)')
ylabel('Pressure (bara)')

%% Display Results
fprintf('Max pressure during the observed time interval = %.3f bara.\n', max(Psol))

%% Functions
    % Main Function
function dUdt = cfunc(t,U_steam)

% Variables
global m_ExtraSteam U0_shell V rho_0 U_steam

Ptube_in = 7;
Ptube_out = 6.5;

T_tube_out = 98;

if t<3.3
m_tube = 707.9;
T_tube_in = 91.7;
else 
m_tube = 1131.2;
T_tube_in = 80.1;
end

% Find the Temp with Internal Energy 
T = fsolve(@TfinderwU,100);

deltaH_tube = XSteam('h_pT', Ptube_out, T_tube_out) - XSteam('h_pT', Ptube_in, T_tube_in);        %  Specific enthalpy difference between the inlet and outlet of the tube side, kJ/kg

if ((m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T))) >= m_ExtraSteam
        m_ExtraSteam_Condensation = m_ExtraSteam;
else
   m_ExtraSteam_Condensation = (m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T)); 
end

m_steam_accumulated = m_ExtraSteam - m_ExtraSteam_Condensation;                 % Steam accumulated in the shell side, kg/sec

steam_mass = rho_0*V + m_steam_accumulated*t;
U_condensate = XSteam('uL_T', T);
CpL = XSteam('CpL_T', T);

Qcv = m_tube*CpL*(T_tube_out-T_tube_in);

% Differential Equation
dUdt = (-Qcv + m_ExtraSteam*U0_shell - m_ExtraSteam_Condensation*U_condensate - U_steam*m_steam_accumulated)/steam_mass;

end

    % Temperature Finder Function with Internal Energy 
function Obj = TfinderwU(T)

global U_steam

uV_T = XSteam('uV_T', T);

Obj = U_steam- uV_T;

end

[Stephen23]

dUdt is empty when the error occurs. Possible causes (that should be adressed anyway):

• Using GLOBAL (with default = empty) instead of parameterizing the function calls: https://www.mathworks.com/help/matlab/math/parameterizing-functions.html. If a variable is used in multiple functions then use nested functions.
• using matrix operations everywhere, a complete lack of array operations: https://www.mathworks.com/help/matlab/matlab_prog/array-vs-matrix-operations.html. This makes it unlikely that the function CFUNC was written to fulfill the requirements given in the ODE45 documentation: "The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector..." source: https://www.mathworks.com/help/matlab/ref/ode45.html#bu00_4l_sep_shared-odefun.

Lets test it right now:

global m_ExtraSteam U0_shell V rho_0 U_steam
T0_shell = 105;                                       % initial shell side temp, degC
P0_shell = 1.2;                                  % initial pressure, bara
U0_shell = XSteam('uV_p', P0_shell);
Plimit  = 3.3;                                        % max pressure, bara
MW      = 0.018015;                             % kg/mol
R       = 0.00008205736;                        % bar.m3/mol.K
V       = 233.21;                                   % m3
rho_0   = XSteam('rhoV_p', P0_shell); 
tb      = 20;                                  % Upper time interval limit
m_ExtraSteam  = 35.08;           % Flowrate of extra steam due to ST trip, kg/sec

cfunc(1,[1,2]) % fails the requirements

Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.

Equation solved at initial point.

fsolve completed because the vector of function values at the initial
point is near zero as measured by the value of the function tolerance,
and the problem appears regular as measured by the gradient.


ans =

     []

function dUdt = cfunc(t,U_steam)
% Variables
global m_ExtraSteam U0_shell V rho_0 U_steam
Ptube_in = 7;
Ptube_out = 6.5;
T_tube_out = 98;
if t<3.3
m_tube = 707.9;
T_tube_in = 91.7;
else 
m_tube = 1131.2;
T_tube_in = 80.1;
end
% Find the Temp with Internal Energy 
T = fsolve(@TfinderwU,100);
deltaH_tube = XSteam('h_pT', Ptube_out, T_tube_out) - XSteam('h_pT', Ptube_in, T_tube_in);        %  Specific enthalpy difference between the inlet and outlet of the tube side, kJ/kg
if ((m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T))) >= m_ExtraSteam
        m_ExtraSteam_Condensation = m_ExtraSteam;
else
   m_ExtraSteam_Condensation = (m_tube*deltaH_tube)/(XSteam('hV_T', T) - XSteam('hL_T', T)); 
end
m_steam_accumulated = m_ExtraSteam - m_ExtraSteam_Condensation;                 % Steam accumulated in the shell side, kg/sec
steam_mass = rho_0*V + m_steam_accumulated*t;
U_condensate = XSteam('uL_T', T);
CpL = XSteam('CpL_T', T);
Qcv = m_tube*CpL*(T_tube_out-T_tube_in);
% Differential Equation
dUdt = (-Qcv + m_ExtraSteam*U0_shell - m_ExtraSteam_Condensation*U_condensate - U_steam*m_steam_accumulated)/steam_mass;
end
    % Temperature Finder Function with Internal Energy 
function Obj = TfinderwU(T)
global U_steam
uV_T = XSteam('uV_T', T);
Obj = U_steam- uV_T;
end

[Talha]

Thanks. Most likely, this happens because of the global variable "U_steam" which is empty at the start.