%direction field logistic equation script
close all
clear all

a = 2; %params for logistic
b = 1; %ditto
f = @(N) (N.*(a-b*N)); %right hand side of DE as a fctn handle
N0_1 = 0.75; %starting values
N0_2 = 4; %ditto
N0_3 = 0.1; %ditto again
N_true = @(t,N0) (a/b)./( 1 + ( (a-b*N0)/(b*N0) )*exp(-a*t)  ); %true soln as fctn handle

h = 0.2; %step size for vector field
h_soln = 0.05; % step size for true solution (should be smaller than h)
[t,N]=meshgrid(0:h:4,0:h:4); % meshgrid domain
tp = [0:h_soln:4]; % time vector for true sol'n

N_true1 = N_true(tp, N0_1); % true soln for N0_1
N_true2 = N_true(tp, N0_2); % true soln for N0_2
N_true3 = N_true(tp, N0_3); % true soln for N0_3

dN = f(N);
dt = ones(size(dN));
dNu = dN./sqrt(dt.^2+dN.^2);
dtu = dt./sqrt(dt.^2+dN.^2);
%dNu = dN; %non-normalized
%dtu = dt; %ditto

figure(1)
quiver(t,N,dtu,dNu,'m'); % plot direction field normalized
hold on;
plot(tp,N_true1,'b-',tp,N_true3,'r-',tp,N_true2,'k-','LineWidth',1.25) %sol'n plots
hold off;

legend('Direction Field','solution for N0<a/b','solution for N0<a/b','solution for N0>a/b')


axis([0 4, 0 4])
xlabel('time: t')
ylabel('Population: N(t)')
title('Stunning Direction Field Plot')