Compute the basin of attraction for a steady state of an ordinary differential equation ODE

I need write a function called "compute_basin" that takes in a steady state and its k value and computes which set of initial conditions in [15,30] approach that steady state. Specifically it will return two values, "a" and "b", which describe the set [a,b] of initial conditions that approach the steady state. The function will be assessed by running
[a,b]=compute_basin(steady_state,k);
for many different values of k and the associated steady states.
Inputs:
a scalar value of a steady state
a scalar value of k (note the steady state should be correct for the k value)
Outputs:
a scalar for the lower bounds of the basin of attraction
a scalar for the upper bounds for the basin of attraction
This is what i have so far:
function [a, b] = compute_basin(Ceq, k)
% Initialize variables
a = 15; % Lower bound for the basin of attraction
b = 30; % Upper bound for the basin of attraction
threshold = 0.001; % Error threshold

% Binary search for the lower bound (a)
left = 15;
right = 30;
while left < right
mid = (left + right) / 2;
tstar = compute_time_to_steady(mid, k, Ceq);
if tstar > 0
right = mid;
else
left = mid + 0.1;
end
end
a = left;

% Binary search for the upper bound (b)
left = 15;
right = 30;
while left < right
mid = (left + right) / 2;
tstar = compute_time_to_steady(mid, k, Ceq);
if tstar > 0
left = mid;
else
right = mid – 0.1;
end
end
b = right;
end
But unfortunalty this is not correct. I have a tip saying that "There are many ways of doing this but a binary search algorithm might be helpful (Wikipedia entry for binary searchLinks to an external site.). Make sure that the [a,b] bounds are within the range of [15,30]. Note that some steady states may be unstable and so their basin of attraction will only be their exact value. In these cases a=b. It helps to solve for the upper and lower bounds separately, i.e. running the code for a and then running it again for b."

I have acces to three functions, the first is "compute_ode" that computes the derivative dC/dt for the differential equation model dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k.
function dcdt = compute_ode(t,c,k)
% write the equation for the differential equation and assign it to output dcdt
dcdt = 0.1 * (c-20) * (23-c) * (c-26) + k;

end
Then i have the second "compute_states" that takes as input a value for k and then returns as output all real-valued steady states for the differential equation.
function vectCeq = compute_states(k)
% Coefficients of the polynomial equation -0.1*c^3 + 6.9*c^2 – 157.8*c + 1196 + k = 0
coeffs = [-0.1, 6.9, -157.8, 1196 + k];

% Find all roots of the polynomial
all_roots = roots(coeffs);

% Filter out complex roots
real_roots = all_roots(imag(all_roots) == 0);

% Return the real-valued steady states
vectCeq = real_roots;
end
And latsly, I have "compute_time_to_steady" that solves the differential equation and returns the earliest time it takes to get within 1% error of the steady state value:
function tstar = compute_time_to_steady(init_cond, k, steady_state)
% Set the options for the ODE solver
options = odeset(‘RelTol’, 1e-6, ‘AbsTol’, 1e-9);

% Define the ODE function as a nested function
function dcdt = ode_function(t, c)
dcdt = compute_ode(t, c, k);
end

% Initialize variables
tspan = [0 100000]; % Time range from 0 to 100,000
threshold = 0.001; % 1% error threshold

% Solve the ODE
[T, C] = ode45(@ode_function, tspan, init_cond, options);

% Iterate over the solution to find the first time it meets the criterion
for i = 1:length(T)
if abs(C(i) – steady_state) / steady_state < threshold
tstar = T(i);
return;
end
end

% If the criterion is not met within the time range, return the last time point
tstar = T(end);
end

function dcdt = compute_ode(t, c, k)
dcdt = 0.1 * (c – 20) * (23 – c) * (c – 26) + k;
end

function steady_states = compute_states(k)
% Define the steady state values based on the value of k
steady_states = [20; 23; 26];
end

I would really appriciate any help on computing the basin of attraction for a steady state!
Thank you!I need write a function called "compute_basin" that takes in a steady state and its k value and computes which set of initial conditions in [15,30] approach that steady state. Specifically it will return two values, "a" and "b", which describe the set [a,b] of initial conditions that approach the steady state. The function will be assessed by running
[a,b]=compute_basin(steady_state,k);
for many different values of k and the associated steady states.
Inputs:
a scalar value of a steady state
a scalar value of k (note the steady state should be correct for the k value)
Outputs:
a scalar for the lower bounds of the basin of attraction
a scalar for the upper bounds for the basin of attraction
This is what i have so far:
function [a, b] = compute_basin(Ceq, k)
% Initialize variables
a = 15; % Lower bound for the basin of attraction
b = 30; % Upper bound for the basin of attraction
threshold = 0.001; % Error threshold