M = [1,0;0,0];
a = 100;

tol = 1.e-6
options_dae = odeset('Mass',M,'AbsTol',tol,'RelTol',tol,'Stats','on');
options_ode = odeset('AbsTol',tol,'RelTol',tol,'Stats','on');

disp('DAE')
[t_dae,y_dae] = ode15s(@(t,y) dae(t,y,a), [0,10*pi],[0;0],options_dae);
disp(''), disp('SSP')
[t_ssf, y_ssf] = ode15s(@(t,y) ssf(t,y), [0,10*pi],[0],options_ode);
disp(''),disp('ODE')
[t_ode, y_ode] = ode15s(@(t,y) ode(t,y,a), [0,10*pi],[0;0],options_ode);


a = 100;
tol = 1.e-2;
err = 1;
disp('DAE')
while err > 1.e-6
  tol = tol/2;
  options_dae = odeset('Mass',M,'AbsTol',tol,'RelTol',tol,'Stats','off');
  solution = ode15s(@(t,y) dae(t,y,a), [0,10*pi],[0;0],options_dae);
  err = norm(solution.y(:,end),inf);
end
tol
err
solution.stats
disp('')

disp('SSF')
tol = 1.e-2;
err = 1;
while err > 1.e-6
  tol = tol/2;
  options_ode = odeset('AbsTol',tol,'RelTol',tol,'Stats','off');
  solution = ode15s(@(t,y) ssf(t,y), [0,10*pi],[0],options_ode);
  err = norm(a*(solution.y(:,end))^2,inf);
end
tol
err
solution.stats
disp('')

disp('ODE')
tol = 1.e-8;
err = 1;
while err > 1.e-6
  tol = tol/2;
  options_ode = odeset('AbsTol',tol,'RelTol',tol,'Stats','off');
  solution = ode15s(@(t,y) ode(t,y,a), [0,10*pi],[0,0],options_ode);
  err = norm(solution.y(:,end),inf);
end
tol
err
solution.stats


function dy = dae(t,y,a)
	dy = [y(2)-a*y(1)^2+cos(t); -y(2)+a*y(1)^2];
end

function dy = ssf(t,y)
	dy = cos(t);
end

function dy = ode(t,y,a)
	dy = [y(2)-a*y(1)^2+cos(t); 2*a*y(1)*(y(2)-a*y(1)^2+cos(t))];
end