%%%%Script for solving the ODE x' = f(t,x) = ax(1-x) + h(1+sin(2pi t))
close all
clear all
%define constants
h = 0.8; %0.8, 0.91,1.1 1.2
a = 4;

%%initial condition(s)
X0 = [0.32:0.2:2];

%initial time
ti = 0;

%final time 
tf = 3;

tspan = [ti tf];

%define right-hand side function of ODE:
fun = @(tt,xx) a.*xx.*(1-xx) - h.*(1+sin(2.*pi.*tt));


%%make a slope field
[tm,xm] = meshgrid(-2:0.05:2);
dx = 1;
dy = fun(tm,xm);
r = ( dx.^2 + dy.^2 ).^0.5;
px = dx./r;
py = dy./r;
figure(10)
quiver(tm,xm,px,py);
xlim([-2 2])
ylim([-2 2])

pause


%obtain solution vectors for the initial conditions x0, on the time range
%ti to tf:
for ii = 1:length(X0)
    x0 = X0(ii);
    [tt,xx] = ode45(fun, tspan, x0);  %xx gives the vector of solution values at the corresponding time points in tt
    
    figure(1)
    hold on
    plot(tt,xx,'o-','LineWidth',3)
    hold off
    xlabel('t')
    ylabel('x(t)')  
    xlim([0 tf])
    ylim([0 2])
end

pause

%%%% Now computationally "sketch" the poincare map p, take a bunch of
%%%% initial conditions, and plot them against the corresponding solutions
%%%% at t = 1;


X00 = [0.0:0.005:1.25];
X1 = zeros(1,length(X00));

%initial time
ti = 0;

%final time 
tf = 1;

tspan = [ti tf];

%define right-hand side function of ODE:
fun = @(tt,xx) a.*xx.*(1-xx) - h.*(1+sin(2.*pi.*tt));


%%make a slope field
[tm,xm] = meshgrid(-2:0.05:2);
dx = 1;
dy = fun(tm,xm);
r = ( dx.^2 + dy.^2 ).^0.5;
px = dx./r;
py = dy./r;
figure(10)
quiver(tm,xm,px,py);
hold on
plot(tt,xx,'LineWidth',3)
hold off
xlim([0 2])
ylim([0 2])
xlabel('t')
ylabel('slopes f(t,x)')




%obtain solution vectors for the initial conditions x0, on the time range
%ti to tf:
for ii = 1:length(X00)
    x0 = X00(ii);
    [tt,xx] = ode15s(fun, tspan, x0);  %xx gives the vector of solution values at the corresponding time points in tt
    X1(ii) = xx(end);
    
    figure(2)
    plot(X00, X1,'.','MarkerSize',10)
    hold on
    plot([0:0.1:3],[0:0.1:3],'-.','LineWidth',2)
    hold off
    xlabel('x_0')
    ylabel('p(x_0) = x(1)')
    xlim([-0.5 2])
    ylim([-0.5 2])
       
end

