Implement solution of differential equation within another differential equation

Hello,

I am working on a project that requires the solution of a differential equation to be used within another differential equation, however when I attempt to create a new file I get the error message:

*Error using feval Undefined function 'acceleration' for input arguments of type 'double'.

Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Acceleration_solve (line 4) [t,U]=ode45('acceleration',tspan,U_0)*

Is there a way to implement the differential equation solution within the new differential equation? I will attach the code.

function U_prime=acceleration(t, U)

%%Solve differential equation and plot (This was defined in another file)
V_0=1.2;                         % Define initial condition for air volume [L]
tspan=[0:0.01:100];              % Define time span for solution [ms]
[t,V]=ode45('vol_air',tspan,V_0); % Solve ODE

idx=find(V<2);                   % Find all indexs where bottle has water
te=t(idx);                       % Find time of air expulsion

plot(t,V)
xlabel('Time [ms]')
ylabel('Volume of air in bottle [L]')

%%Variable definitions
Pa=101325;                       % Atmospheric pressure [Pa]
P0_psi=80;                       % Internal bottle pressure in PSI initial [PSI]
P0=(P0_psi*6894.76)+Pa;          % Convert PSI to Pascal add 1 atm for gaugue pressure [Pa]
V_bottle=2;                      % Volume of bottle [L]
V_water=0.8;                     % Volume of water [L]
V_0=V_bottle-V_water;            % Initial volume of air [L]
rho_w=1000;                      % Density of water [kg/m^3]
gamma=1.4;                       % Specific heat ratio for air (Cp/Ct)
Dn=0.022;                        % Nozzle diameter in metres [m]
An=(pi*(Dn^2))/4;                % Nozzle exit area [m^2]
P=P0.*(V_0./V).^gamma;           % Pressure at time t, P(t) [Pa]
R_gas=287;                       % Gas constant for air     
T0=P0*(273.15/Pa);               % Temperature at t=0, P0 initial PSI pressure [K]
T=P*(T0/P0);                     % Temperature at time t, T(t) [K]
g=9.81;                          % Gravitational acceleration on earth [ms^-2]
Cd=0.2;                          % Coefficient of drag

%%M(t)
M_dry=0.1;                       % Dry rocket mass
M_air=P.*V./((R_gas.*T0).*1000); % Mass of air at time t Mair(t)
% M_air_true=find(M_air>0);      % Check air mass never becomes negative  
M_t=M_dry+M_air+(rho_w.*((2/1000)-(V./1000))); % Total mass at time t, M(t)
M_allowed=M_t(idx)               % Mass of rocket under allowed conditions, i.e. V<=2. 

%%ODE for acceleration
U_prime(1,1)=(((2.*An)./M_t).*(P-Pa))-g-(0.5.*U.^2.*S.*Cd./M_t)

I have then run the ODE45 solver in a different file with the following code:

%%Solve differential equation and plot

U_0=0;
tspan=[0:0.01:100];
[t,U]=ode45('acceleration',tspan,U_0)

%%Plot
plot(t,U)
xlabel('Time [s]')
ylabel('Acceleration')