function doall=sea_ice_switch
clear all; close all;
% Eli Tziperman: example for course in nonlinear dynamics and
% chaos:  Solve simple sea ice switch equation.

% define parameters:
% ------------------
% Initial values:
Temp=10.0; V=0.5; x(1:2,1)=[V,Temp]'; %x(1,1)=Temp; x(2,1)=V;
% time step: needs to be small enough for a given problem:
delta_t=0.1;
% Print (and save) every nprint steps:
nprint=1;
% Integration time:
T=500.0;

% Number of steps:
N=floor(T/delta_t);
NP=floor(N/nprint);


% Integrate:
for n=1:1:NP
    xp=x(:,n)';
    for j=1:1:nprint
      x_tilde=xp+delta_t*rhs(xp);
      xp=xp+0.5*( rhs(xp) + rhs(x_tilde) )*delta_t;
    end
    x(:,n+1)=xp';
    time(n)=n*delta_t;
end

% Plot simple time series:
plot(time(1:n),x(1,1:n))
xlabel('time')
ylabel('ice volume')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff=rhs(x)
% return the rhs of the dynamical system:

V=x(1);
Temp=x(2); 

gamma=-0.1; Td=-2.; 
a1=0.005; a2=0.0025; a3=3.; a4=54.; a5=25.0; a6=64.;

% ice volume equation:
ff(1)=a1+a2*Temp;
% temperature equation:
ff(2)=(-a3*Temp+a4-a5*Theta(Temp-Td)-a6*V);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tt=Theta(X)
% smoothed Heaviside function, one for X<0, zero for X>0:
tt=0.5*(1.0-tanh(X/0.1));

