r = 2;
deltam = 1;
deltap = 0.5;
gammam = 1;
gammap = 1;
alpha = 10;
beta = 1;
a = (alpha*beta)/(deltam*deltap);
k = 0.1;
Tf = 200;
y0 = [0;0;0;1];%y(end,:);
allts=[];
allys=[];
Tstart=0;
disc = a*(a-4*k);
root1 = 1/2*(a-2*k+sqrt(disc));
root2 = 1/2*(a-2*k-sqrt(disc));
% Runs through this 4 times, every 200 seconds 
% changing deltap to give "switch" behavior
for deltap = [0.5 2 0.5 2]
odefun = @(t,y) CircadianRhythm(t,y,alpha,beta,gammam,gammap,deltam,deltap,r,k);
[t,y] = ode45(odefun, [0 Tf], y0);
allts=[allts;t+Tstart];
Tstart=Tstart+Tf;
allys=[allys;y];
y0=y(end,:);
end
plot(allts,allys(:,4))

function dydt = CircadianRhythm(t,y,alpha,beta,gamma_m,gamma_p,delta_m,delta_p,r,k)
    mn = y(1);
    mc = y(2);
    pc = y(3);
    pn = y(4);
    dydt = [alpha*(pn/(k+pn))^r+1e-2 - gamma_m*mn;...
        gamma_m*mn - delta_m*mc; ...
        beta*mc-gamma_p*pc; ...
        gamma_p*pc - delta_p*pn];
end
