Solution Plot of the Constant Harvesting Function with Directional Field

I'm trying to plot solutions along with the directional field as in the image given below, of the constant harvesting model given by

dP/dt = r*P(1 - (P/M)) - h , P(0) =p0 , where r = 0.8, M = 8000 and h = 0 in this particular case.

I attempted using the following code, which does not result in something as in the image given below. I would appreciate your help with making the graph looks like as in the image above.

% Constant Harvesting Model

% Parameters
r = 0.8;    % Intrinsic growth rate
K = 8000;   % Carrying capacity

% Define the system of differential equations
dN_dt = @(N) r * N .* (1 - N/K);

% Define the range of values for N and t
N_range = 0:2000:14000;
t_range = 0:5:20;

% Create a grid of N and t values
[N, t] = meshgrid(N_range, t_range);

% Compute the rate of change of N at each point on the grid
dN = dN_dt(N);

% Create a figure
figure;
hold on;

% Plot the directional field
quiver(t, N, ones(size(dN)), dN, 'k');

% Set plot properties
xlabel('t');
ylabel('N');
title('Constant Harvesting Model');
legend('Directional Field');
grid on;

% Invert the t-axis
set(gca, 'XDir', 'reverse');

% Display the figure
hold off;

Answer 1

[Sam Chak]

Hi @Gith

I'm unsure, but I think that it was the quiver scaling issue when the carrying capacity K is too large (on the Y-axis), compared to the time scale (on the X-axis). One workaround is to rescale on the Y-axis, and in your case, it is related to the value K.

r = 0.8;
K = 8000/1000;
[T, N] = meshgrid(0:20/14:20, 0:K/14:K);
S = r*N.*(1 - N/K);
L = sqrt(1 + S.^2);
U = 1./L;
V = S./L;
quiver(T, N, U, V, 0.3)
axis square

hold on

% differential equation
f = @(T, N) r*N.*(1 - N/K);

tspan  = 0:0.01:20;
init   = 10/1000;
[T, N] = ode45(f, tspan, init);
plot(T, N, 'r', 'linewidth', 1.5)

hold off
grid on
xlabel({'$t$'}, 'Interpreter', 'latex')
ylabel({'$N(t), \quad (\times 1000)$'}, 'Interpreter', 'latex')