% Logistic Model
% dP/dt=k*P*(1-P/N);


% Clear all variables

clear all

% Slope field

[t,p]=meshgrid([0:0.5:5], [0:0.5:6]);
m=f(t,p);
L=sqrt(1+m.^2);
quiver(t,p,1./L,m./L,'.r');
axis([0 5 0 6])
hold on;

% Numerical solution using ode45

t0=0;h=0.1;tfinal=5;p0=0.5;
[x,y]=ode45(@(t,p)f(t,p),[t0:h:tfinal],p0);
plot(x,y,'-','LineWidth',2);
hold on

%Numerical solution by using ode1=Euler's method

t0=0;h=0.5;tfinal=5;p0=0.5;
zout=ode1(@(t,p)f(t,p),t0,h,tfinal,p0);
plot([t0:h:tfinal],zout,'b-o','LineWidth',1)
hold off


% Code realizing Euler's method

function zout=ode1(F,t0,h,tfinal,y0)
z=y0;
zout=z;
for t=t0:h:tfinal-h
    s=F(t,z);
    z=z+h*s;
    zout=[zout;z];
end
end

% The function on the right-hand-side of ODE

function y=f(t,p)
y=2*p.*(1-p/4);
end



