function dy = foxrabbit(t,y);
dy = zeros(2,1);
global a b c d w w1 w2
% y(1) = # of rabbits
% y(2) = # of foxes
dy(1) =  a*y(1) - b*y(2)*y(1);
dy(2) = -d*y(2) + c*y(1)*y(2);
% add climate change oscillation cos(w*t)
dy(1) = a*y(1)*(1+w1*cos(w*t))- b*y(2)*y(1);
dy(2) = -d*y(2)*(1+w2*cos(w*t)) + c*y(1)*y(2);

