Population Growth Model Analysis

Below I tried to solve Population Growth problem, I provide you the mathematical solution. Also, I create MATLAB code, however i feel like that i repeat this code to each part of this code as you can see. How could I improve this code?

Assume a continuous growth model where the function that gives the rate of change is defined as . Provide the function , which represents the population size.
If the growth of a population follows the logistic law but there is also a term for population reduction (e.g., from a predator while the members of the population are the prey), the equation that expresses the population size is of the form
where the term is given here to be and it represents a term for population reduction. The constants are positive. Normalize the equation by setting and .Find the steady solutions of the problem (real and positive) and then implement a program in MATLAB to numerically solve the dimensionless form of the equation. Make the graphical representation of the numerical solution for , with initial conditions , and respectively. What do you observe?

So we have

We are going to find the function :
To find , we integrate :

where is the integration constant, which can be determined if an initial condition is given, such as .
Now we are going to normalize the equation:
Let and . Then:

We substitute into the original equation:

Then we are going to calculate the Stationary Points:
Set :
]

For :

% Define the normalized logistic equation with reduction term
% K/a = 35, ar/b = 0.4
Ka = 35;
ar_b = 0.4;

% Define the function for the normalized logistic equation
f = @(tau, N) ar_b * N * (1 – N/Ka) – N / (1 + N);

% Time span for the simulation
tspan = [0 50]; % Adjust if needed to observe long-term behavior

% Initial conditions
N0 = [40, 25, 2, 0.05];

% Prepare the figure for plotting
figure(1);
hold on;

% Colors for different trajectories
colors = [‘r’, ‘g’, ‘b’, ‘k’]; % Red, Green, Blue, Black

% Solve the equation for each initial condition and plot
for i = 1:length(N0)
[T, N] = ode45(f, tspan, N0(i));
plot(T, N, ‘Color’, colors(i), ‘LineWidth’, 2);
legendInfo{i} = [‘N_0 = ‘ num2str(N0(i))]; % Create legend entry
end

% Add graph details
title(‘Numerical solutions of the normalized logistic equation’);
xlabel(‘Normalized time τ’);
ylabel(‘Normalized population size N’);
legend(legendInfo, ‘Location’, ‘northeastoutside’);
grid on;
hold off;

% Equilibrium analysis
syms N
eqn = ar_b * N * (1 – N / Ka) – N / (1 + N) == 0;
equilibria = double(solve(eqn, N));

% Stability analysis
f = @(N) ar_b * N * (1 – N / Ka) – N / (1 + N);
J = diff(f(N), N);
eigenvalues = double(subs(J, N, equilibria));

disp(‘Equilibrium Points and their Stability:’);
for i = 1:length(equilibria)
if real(eigenvalues(i)) < 0
stability = ‘Stable’;
else
stability = ‘Unstable’;
end
fprintf(‘N = %f: Eigenvalue = %f, %sn’, equilibria(i), eigenvalues(i), stability);
end

% Bifurcation analysis setup
ar_b_values = linspace(0.1, 1, 100);
bifurcation_results = zeros(length(ar_b_values), 10); % Adjust the columns based on expected number of roots

% Compute equilibrium points for varying ar/b
for i = 1:length(ar_b_values)
eqn = ar_b_values(i) * N * (1 – N / Ka) – N / (1 + N) == 0;
solutions = double(solve(eqn, N));
bifurcation_results(i, 1:length(solutions)) = solutions’;
end

% Plot bifurcation diagram
figure;
plot(ar_b_values, bifurcation_results, ‘b*’);
title(‘Bifurcation Diagram’);
xlabel(‘ar/b’);
ylabel(‘Equilibrium Points’);

% Set up the ODE function for phase plane analysis
f = @(N) ar_b * 0.4 * N .* (1 – N / Ka) – N ./ (1 + N);

% Create a grid of N values and compute their time derivatives
N_values = linspace(0, 50, 400);
dN_dtau = f(N_values);

% Plot phase plane
figure;
plot(N_values, dN_dtau);
title(‘Phase Plane Analysis’);
xlabel(‘Population Size N’);
ylabel(‘dN/dτ’);
grid on;
hold on;

% Mark zero crossings (equilibrium points)
zero_crossings = N_values(abs(dN_dtau) < 1e-3);
for i = 1:length(zero_crossings)
plot(zero_crossings(i), 0, ‘ro’);
end
hold off;Below I tried to solve Population Growth problem, I provide you the mathematical solution. Also, I create MATLAB code, however i feel like that i repeat this code to each part of this code as you can see. How could I improve this code?