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); [sol] = ode15s (@rob, [0, 100], [1; 0; 0], opt); sol.te sol.ye sol.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