function dy=model_fun(t,y)
dy = zeros(2,1);    % a column vector
% time step = one year
% constant definition
tau_atm=2.82;      % years
tau_ocean=100.3;   % years
alfa=3.23;         % coupled strength
fa=0.367;          % climate sensitivity for atmosphere [K/W/m^2]
fo=0.0689;         % climate sensitivity for ocean [K/W/m^2]
RF=-20.0;
dy(2)=
dy(1)=dF/Cm-y(1)/(Cm*lambda)-Cd/Cm*dy(2);
dy(2)=1/tau_ocean*(-y(2)+fo*RF);