r = 5;
v = 1e-2;
k = 1e-2;
alphaBeta = 10*(v^2*k*4*5^5);
alpha = 1;
Tf = 3600*24;
y0 = y(end,:);%[0;0;0;0];%y(end,:);
Beta = alphaBeta/alpha;
odefun = @(t,y) CircadianRhythm(t,y,alpha,Beta,v,r,k);
[t,y] = ode45(odefun, [0 Tf], y0);


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