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 1. I would appreciate your help with making the graph looks like as in the image 1.

Image 1(Required graph): drive[DOT]google[DOT]com/file/d/1WPREaQrGLi_cGrXevtN5ajMq1QR4acGI/view?usp=sharing

Image 2(The result I got with my code): drive[DOT]google[DOT]com/file/d/1QkDFIsTO0UzZfRCpBhfZtiPnpqlWPxL0/view?usp=sharing

I'm using MATLAB online free version

Unfortunately, as an AI language model, I am not able to access external links or images. However, I can suggest some modifications to your code that might help you achieve the desired graph. Here's what you can try:

1. Define the range of P and t values:
```
P=linspace(0,8000,50);
t=linspace(0,50,50);
```

2. Define the function:
```
r=0.8;
M=8000;
h=0;
f=@(t,P) r*P.*(1-P./M)-h;
```

3. Create a meshgrid to plot the directional field:
```
[T,P]=meshgrid(t,P);
dP_dt=f(T,P);
quiver(T,P,ones(size(dP_dt)),dP_dt);
```

4. Solve the ODE using the initial condition and plot the result:
```
p0=5000; % example initial condition
[t,P]=ode45(f,t,p0);
plot(t,P);
```

Put everything together and you should have a graph similar to the one you want. Note that you might need to adjust the number of points in the linspace and the size of the arrows in the quiver plot to get a better fit.

I'm sorry, but as a text-based AI, I cannot view or access external images or files. However, I can guide you on how to create a directional field plot and solutions graph for the constant harvesting model using MATLAB. Here's a step-by-step approach:

1. Import the required constants and initial conditions:
```matlab
r = 0.8;
M = 8000;
h = 0;
p0 = <initial_value>; % specify the initial population value
```

2. Define the differential equation:
```matlab
dPdt = @(t, P) r * P * (1 - (P / M)) - h;
```

3. Create a meshgrid for the P and t axes:
```matlab
tspan = [0 500]; % specify the time range for the solutions
P = linspace(0, M, 30); % specify the population range
[T, P] = meshgrid(tspan, P);
```

4. Use quiver function to create the directional field plot:
```matlab
dP = dPdt(T, P);
quiver(T, P, ones(size(dP)), dP, 'k');
```

5. Use ode45 function to numerically solve the differential equation:
```matlab
[t, y] = ode45(dPdt, tspan, p0);
```

6. Plot the solutions:
```matlab
plot(t, y, 'r', 'LineWidth', 2);
xlabel('Time');
ylabel('Population');
legend('Directional Field', 'Solutions');
title('Constant Harvesting Model');
```

Make sure to replace `<initial_value>` with the desired initial population value. You can adjust the parameters and settings as per your requirements.

Please note that without being able to see the intended graph or your current code, I can only provide a general solution based on the information you provided.