fsolve vs. running a model to equilibrium

Hey all,

So I've been asking a lot of questions about a set of models I'm trying to find the equilibrium of and I've gotten a lot of great help on here (thank you!) but I'm still having trouble. I'm using fsolve to find the equilibrium values of a series of four differential equations (with the code I've copied below) and using this command in the command window:

>> x = fsolve(@odefun_dynamic, [5, 0.000001, 0.1420, 1e-7], options)

And I'm getting very different equilibrium values from when I just run the model to equilibrium and see what values it ends at. I think I could be hitting other local minima instead of the global, but I've been varying the x0 values for the fsolve command and I still don't get an equilibrium consistent to the one I see when running the model (which is the same set of differential equations). Are there other options I should set besides TolFun? Should I try using a different function than fsolve? Any help would be much appreciated! Thank you!!

function x_prime = odefun_dynamic(x)

options = optimset('TolFun',1e-14);

% Flow rate parameter
F = .01; % flow rate of medium

%set parameters
S = 1e-7; %initial substrate concentration
p_max = 3.38e-6; %max intake rate of nutrient (P) by algae
K_p = 1.29e-8; %half saturation constant of nutrient intake by algae
u_Amax = 1; %apparent max growth rate
Q_Amin = (2.5E-7); %subsistence quote for algae
q_max = 1; %max uptake rate of algae by Daphnia
K_q = 0.164; %half saturation constsant of uptake of algae by Daphnia
gamma = 0.5; %fraction of biomass of algae assimilated by Daphnia
m_d = 0.02; %mortality of Daphnia

% Get state variables from x
N_A = x(1);
Q_A = x(2);
N_d = x(3);
R = x(4);

%functions

M_A = F * N_A;

p = (p_max * R) / (K_p + R);

I_R = N_A * p;

r_A = (u_Amax * N_A) * ((Q_A - Q_Amin) / (Q_A));

M_d = m_d * N_d;

I_A = (N_d * q_max * N_A) / (K_q + N_A);

r_d = gamma * I_A;

% differential equations

N_A_prime = r_A - M_A - I_A;

Q_A_prime = p - ((Q_A - Q_Amin) * u_Amax);

N_d_prime = r_d - M_d;

R_prime = (F * S) - (F * R) - I_R;

% Return derivatives in a single column vector
x_prime = [N_A_prime; Q_A_prime; N_d_prime; R_prime];

end