The ode45 function only outputting constant value

I am currently trying to create a program which accepts user inputs for parameters, passes those to a function, and then passes that function to the ode45 solver and plots the results. However, the ode solver seems to only give a constant output equal to the specified initial conditions. I was wondering whether there were some kind of error in my script that was leading to this result.

%% Units: g, mol, dm^3, K, Kpa, s
clear
clc
%Specifying operating conditions
    R = 8.314472; %Ideal Gas Constant
    m_a = 60.052; % Molar mass acetic acid
    m_b = 74.1; % Molar mass Butanol
    phi = 0.2; %Porosity of Catalyst 
    rho_c = 19300; %Density of Catalyst
    P_0 = input('Please specify an operating pressure in Kpa.');
        while P_0 <= 0
            P_0 = input('Please specify an operating pressure greater than zero.');
        end
    T = input('Please specify an operating temperature in K.');
        while T <= 0
            T = input('Please specify an operating temperature greater than zero.');
        end
    Dp = input('Please specify a pipe diameter in dm.');
        while Dp <= 0
            Dp = input('Please specify a pipe diameter greater than zero.');
        end
    F_t_0 = input('Please specify incoming molar flow rate in mol/s.');
        while F_t_0 <= 0
            F_t_0 = input('Please specify an incoming molar flow rate greater than zero.');
        end
    x_b = input('Please specify mole fraction of Butanol in feed.');
        while x_b <= 0 || x_b > 1
            x_b = input('Please specify a mole fraction between zero and one.');
        end
    x_a = 1 - x_b; %Mole fraction of acetic acid
%Calculating parameters for pressure drop equation
    Ap = pi*(Dp/2)^2; %Cross Sectional Area of Pipe
    Cto = P_0/(R*T); %Initial concentration
    rho_a = x_a*Cto*m_a; %Density of acetic acid
    rho_b = x_b*Cto*m_b; %Density of butanol
    rho = x_a*rho_a + x_b*rho_b; %Approximate total density
    mu_a = 0.001155; %Viscosity of acetic acid
    mu_b = 0.002593; %Viscosity of butanol
    mu = x_a*mu_a + x_b*mu_b; %Approximate total viscosity
    m_flow = x_a*F_t_0*m_a + x_b*F_t_0*m_b; %Incoming mass flow rate
    G = m_flow/Ap; %Mass velocity 
    beta_1 = (G*(1-phi)/(rho*Dp*phi^3));
    beta_2 = ((150*(1-phi)*mu)/Dp) + 1.75*G;
    beta = beta_1*beta_2;
    rho_bulk_c = rho_c*(1 - phi); %Bulk density of Catalyst
    alpha = ((2*beta)/(Ap*rho_bulk_c*P_0));
%Calculating parameters for rate law
    k_f = 0.0090953*exp(-45.49/(R*T));
    k_r = 0.008643*exp(-23.90/(R*T));
%ODE solver
w_span = [0:1:100]; %In grams
F_a_0 = F_t_0*x_a;
F_b_0 = F_t_0*x_b;
[W,F] = ode45(@(w,y) f(w,y,Cto,k_f,k_r,F_t_0,alpha), w_span,[F_a_0; F_b_0; 0; 0; 1;]);
plot(W,F(:,1),W,F(:,2),W,F(:,3),W,F(:,4),W,F(:,5))
title('Component Flow Rates and Pressure Drop with Respect to Catalyst Mass');
xlabel('Mass W (g)');
legend('Fa','Fb','Fc','Fd','p')
function dydw = f(w,y,Cto,k_f,k_r,F_t_0,alpha)
   dydw(1) = (Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4))^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4)); %Mole Balance for Species A
    dydw(2) = (Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4))^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4)); %Mole Balance for Species B
    dydw(3) = -1*(Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4)^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4))); %Mole Balance for Species C
    dydw(4) = -1*(Cto^2*y(5)^2/(y(1)+y(2)+y(3)+y(4))^2)*(k_f*y(1)*y(2) - k_r*y(3)*y(4)); %Mole Balance for Species D
    dydw(5) = (-1*alpha/2*y(5))*((y(1)+y(2)+y(3)+y(4))/F_t_0); %Pressure Drop Equation
    dydw = dydw(:);
end   

    

[Walter Roberson]

Please indicate some suitable inputs to test with.

[Braydon Lantz]

For the first four inputs, any positive integer would be acceptable, and for the last, any decimal value between 0 and 1. However, in the tests that I've done, I obtain a constant value regardless of which inputs are used. The intitial conditions used by the ODE solver are defined by the last two inputs, F_t_0 and x_b.

Answer 1

[Walter Roberson]

exp(-45.49/R*T) is exp(-1500-ish) which underflows to 0. Your input temperature would have to be below about 136 K in order for that expression to not underflow to 0.

Perhaps those should be exp(-45.49/(R*T)) ?

[Braydon Lantz]

That was indeed a mistake on my part, and I have made the correction you suggested. However, I am still getting a constant value from the ODE solver. I have edited the script in my initial post to reflect this change.

[Walter Roberson]

What are some realistic inputs? With the random values I am inputting, I am seeing slight curves.

[Braydon Lantz]

That was the problem! I wasn't considering the physical properties of the material when I defined the range of values that are accepted as inputs, so certain numbers were giving anomalous results. Thank you!