Using an ODE to solve for an equation

Hello all,

I am looking to solve a differential equation using ODE45. I have gotten pretty far in my script and having an issues that I imagine is really easy to solve. I have attached the equation I am trying to solve for. Labeled pC_eq. I was able to use the ode function to solve for dVc/dt and rho_c/dt ( I have also attached my function I used to calculate that. It is is called test_fnc.m. I am having a hard time calculating dpc/dt (which is the change in chamber pressure over time). dpc/dt is the last line in my main script below. For my assignment I need to calculate the equlibrium chamber pressure and the chamber pressure. I calculated the equlibrium chamber pressure, I labeled it as Pc_eq in my script)

What I am suppose to see is dpc/dt just slightly of a lesser value of Pc_eq, but that is not happening for me in my last equation.

clear all
%% Step 1
% ========================================================================= %
% Basing on the motor and propellant data indicated in Table 1 and on the 
% burning surface evolution shown in Figs. 5 and 6 and assuming a constant
% chambre temperature profile compute the following curves without taking
% into account the non-ideal parameters.
%% Step 1.1 
%Find: The motor chamber pressure and the vaccume thrust as a
% function of  both time and web assuming quasi-steady state (equilibrium)
% User inputs - First stage P80 SRM
% ========================================================================= %
mp = 88000; % Propellant Mass [kg]
ms = 7330;  % Structural mass [kg]
l = 10.6;   % length [m]
d = 3.0;    % Diameter [m]
f = 3015;   % Max thrust(vaccum) [kN]
bt = 110;   % Buring time [sec]
isp = 280;  % Specifi Impulse (vaccuum) [sec]


% User inputs - Propellant Ballistic Properties 
% ========================================================================= %
a_0 = 1.847e-05; % Temperature coefficient @ 300 K [m/s * Pae-0.4]
n = 0.4;       % Combustion index
tau = 0.0015;  % Temperature sensitivity [k^-1]
rho_p = 1790;    % Density [kg/m^3]
rho_c = 1;        % Initial chamber density [kg

% User inputs - Propellant thermochemocal properties 
% ========================================================================= %
T_F = 3550;   % Flame temperature [K]
M = 29;       % Molecular mass [kg/kmole]
gamma = 1.13; % Specific heat ratio

% User inputs - Motor geometrical properties 
% ========================================================================= %
d_throat = 0.496; % Throat diameter [m]
e = 16;           % expansion ratio
V_c = 8.6;        % Initial chamber volume [m^3]
v_frac = 0.85;    % volumetric loading fraction (V_c/(V_c+V_p)

% User inputs - Pressurizing Gas Properties 
% ========================================================================= %

p_exit = 1.3e+05;     % [Pa]
T_initial = 300;     % [K]

% constants 
% ========================================================================= %
R = 8314.5;       % gas constant [J/kmol)
% Calculations 
% Find: The motor chamber pressure and the vaccum thrust as a
% ========================================================================= %
a_t = pi* (d_throat/2)^2; % Thrat area [m^2]
cap_gamma = sqrt (gamma * (2/(gamma + 1))^((gamma+1)/(gamma-1))); % capital gamma
c_star = (1/cap_gamma)*sqrt((R*T_F)/M);
r_i = d_throat/2; % Need to come back to this value
%L = 10.6; % Need to come back to this value
%beta = a_0* (rho*a_0*c_star*(1/a_t))^(n/(1-n));    %coeff of constants 
%t = 0:0.5:bt;    % Time vector
%aa = r_i^((2*n-1)/(n-1)); % first part of y eq
%bb = ((2*n-1)/(n-1))* (2*pi*L)^(n/(1-n))*beta*t; % second part of y equation
%y = (aa + bb).^((n-1)/(2*n-1)) - r_i;  %burnt web as a function of time
%s_b = 2*pi*(r_i+y)*L;                   %burnt suface area as a function of y
%K = s_b/a_t;   % Klemmung
%Pc_eq = (a_0 *rho*c_star*K).^(1/(1-n));

tow = (V_c/a_t) * (1/(cap_gamma*c_star));  % characteristic time scale

y_0 = 0;
y_f = 1.075;
nn = 100;   %sampling rate
%y_result = cell2mat(arrayfun(@(e) linspace(y_0, y_f, nn), y_f', 'UniformOutput', false));
%y_result = y_result';
time_interval = [0,bt];
y0 = 0;
[t_result,y_result] = ode45(@(t,y) diff_fnc(t,y,a_0,n,d_throat, rho_p, c_star), time_interval, y0);
[~,Pc_eq] = chamberPressure(t_result,y_result,a_0,n,d_throat, rho_p, c_star);

% Calculating the Force

x = 1:0.0001:7;   % defining a reasonable range for Mach number @ the exit
num =  (gamma + 1) /2;
exponent = (gamma+1)/(2*(gamma-1))
dem = (gamma -1)/2;
y = x.*(num./(1+dem*x.^2).^exponent)-0.02;  %simplifed equation for my expansion ratio to M_exit

index = zerocross(y)
y(index)
plot(x,y)
hold on
plot(x(index),y(index),'o')
M_e = x(index);   % Solving for Mach at the exit
ex_1 = (gamma-1)/gamma ;  % I was having the hardest issue with the () so just wrote the exponent to C_f seperately 
c_f = cap_gamma.*sqrt( 2*gamma./(gamma-1) * (1 - (p_exit./Pc_eq).^ex_1)         ) + (p_exit*e)./Pc_eq;

F_1 = c_f .* a_t .* Pc_eq;

%% plots for motor chamber pressure asa function of both time and web
%%% assuming quassi-steady state (equilibrium)
subplot(1,2,1)
plot(t_result,Pc_eq)
xlabel('time [sec]')
ylabel('Pc_{eq} [Pa]')
title ('(Pc)_{eq} as a Function of Time') 
subplot(1,2,2)
plot(y_result,Pc_eq)
xlabel('Burned Web Thickness [m]')
ylabel('Pc_{eq} [Pa]')
title ('(Pc)_{eq} as a Function of Burned Web Thickness')

%% plots for force as a function of both time and web
%%% assuming quassi-steady state (equilibrium)
figure(2)
subplot(1,2,1)
plot(t_result,F_1)
xlabel('time [sec]')
ylabel('Force_{eq} [N]')
title ('Force_{eq} as a Function of Time') 
subplot(1,2,2)
plot(y_result,F_1)
xlabel('Burned Web Thickness [m]')
ylabel('Force_{eq{ [N]')
title ('Force_{eq} as a Function of Burned Web Thickness')


%% Step 1.2/1.4
% Find The motor chamber pressure and the vacuum thrust as a function of
% both time and web integrating the ODE representing the overall mass
% balance.
% The variation of the mass of the gases in the chamber over time,
% seperating the role of density variation and that of chamber volume
% variation
%the role of density variation is negligible all the time

% Repreat the calculation of point 2. Assuming a propellant initial
% temperature of 315 (+15 K wrt to baseline) and of 285 K (-15 K wrt
% baseline). Highlight and explain the difference in terms of pressure and
% thrust curves and their integrals over time


for T_ref = [285;300;315]  % 3 different reference temperatures [K]
    a = a_0*e.^(tow*(T_initial - T_ref)); 
    %Y0 = [8.6; rho_c; 8.6; rho_c; 8.6; rho_c];
   %[t_result, Y_result] = ode45(@(t,Y) diff_fnc_mb(t,Y,a,n,a_t, rho_p, c_star, T_ref,R, t_result, y_result), time_interval, Y0);

end
 Y0 = [8.6; rho_c]; % initial condition of Vc and rho_c
[t_result_285, Y_result_285] = ode45(@(t,Y) diff_fnc_mb(t,Y,a(1),n,a_t, rho_p, c_star, T_ref(1),R, t_result, y_result), time_interval, Y0);
[t_result_300, Y_result_300] = ode45(@(t,Y) diff_fnc_mb(t,Y,a(2),n,a_t, rho_p, c_star, T_ref(2),R, t_result, y_result), time_interval, Y0);
[t_result_315, Y_result_315] = ode45(@(t,Y) diff_fnc_mb(t,Y,a(3),n,a_t, rho_p, c_star, T_ref(3),R, t_result, y_result), time_interval, Y0);


[dpcdt_300, ~] = ode45(@(t,Y) test_fnc(t,Y,a_0,n,a_t, rho_p, c_star, T_initial,R, t_result, y_result,cap_gamma, V_c), time_interval, Y0);

[darova]

It's too complicated. Hard to check where the mistake is