function [t, y, res] = singpertex1( mu, u0 )

% 
% A singular perturbation problem
%
% Leave the 'MassSingular' property at its default 'maybe' to test the
% automatic detection of a DAE.

    function dy = f( t, y )
        dy = [ -y(1)*y(2)
               1 - y(1)^2 - y(2)^2 ];
    end

M = [ 1, 0
      0, mu ];

options = odeset('Mass',M);

sol = ode15s(@f,[0 5],u0,options);

t = RefineMesh( sol.x, 5 );
[y, dy ] = deval( sol, t );
np = length( t );
res = zeros(size(dy));
for i=1:np,
    res(:,i) = dy(:,i) - f(t(i), y(:,i));
end

end