function ode15s_event_test()

fun = @ (t, y) [-0.04*y(1) + 1e4*y(2).*y(3);
                 0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2;
                 y(1) + y(2) + y(3) - 1];

y0 = [1; 0; 0];
tspan = [0, 4*logspace(-6, 6)];
M = [1, 0, 0; 0, 1, 0; 0, 0, 0];

opt = odeset ("Events", @feve, "Mass", @massdensefunstate, ...
              "MStateDependence", "none");

%sol = ode15s (@rob, [0, 100], [1; 0; 0], opt);
[t, y, te, ye, ie] = ode15s (@rob, [0, 100], [1; 0; 0], opt);

te
ye
ie

y(:,2) = 1e4 * y(:,2);
figure (2);
semilogx (t, y, "o");
xlabel ("time");
ylabel ("species concentration");
title ("Robertson DAE problem with a Conservation Law");
legend ("y1", "y2", "y3");


function res = rob (t, y)
 res = [-0.04*y(1) + 1e4*y(2).*y(3);
         0.04*y(1) - 1e4*y(2).*y(3) - 3e7*y(2).^2;
         y(1) + y(2) + y(3) - 1];
end

function masrob = massdensefunstate (t, y)
 masrob = [1, 0, 0; 0, 1, 0; 0, 0, 0];
end

function [val, isterminal, direction] = feve (t, y, ~)
 isterminal = [0, 1];
 if (t < 1e1)
   val = [-1, -2];
 else
   val = [1, 3];
 end

 direction = [1, 0];
end

end
