% THIS CODE REQUIRES THE STATISTICS/MACHINE LEARNING TOOLBOX TO BE INSTALLED

% here, we integrate Eq. (1) using a simple explicit Euler method 
% where we add random noise at each timestep
% we're also going to change the parameter b linearly as time progresses

% again, functions defined at the end

% DEFINE PARAMETERS: this is what you might need to change
% define constants r and a, and the start and end values of b
r = 1;
a = 8; 
b_start = 14; 
b_end = 15;

% set up start and end times of the integration
tstart = 0;
tend = 250;

% timestep (recommend not touching this)
dt = 0.05;

% amplitude of random noise 
amplitude = 0.5;

% INTEGRATE AND PLOT  
t = linspace(tstart,tend,floor((tend-tstart)/dt));
b_ramp = linspace(b_start,b_end,length(t));

N = zeros(length(t),1);
N(1) = stable_fixed_point(r,a,b_ramp(1));
 
for i=2:length(t)
     N(i) = N(i-1) + model(N(i-1),r,a,b_ramp(i))*dt + normrnd(0,amplitude*dt);
end

plot(t,N)
xlabel('Time')
ylabel('N')

% DEFINE FUNCTIONS
% function to calculate the stable fixed point
function fp = stable_fixed_point(r,a,b)
    fp = a/2 + sqrt(a^2/4-b); 
end

% defining the model: i.e. f(N,r,a,b)
function ndot = model(N,r,a,b)
    ndot = r*N*(-N^2+a*N-b);
end
