Solving a DAE system over a longer time span using ode15s

I'm working with the example code on how to set up a DAE problem and solve it with either ode15i or ode15s as given in the documentation: http://se.mathworks.com/help/symbolic/set-up-your-dae-problem.html?refresh=true

However, I'm interested in solving the set of equations for a longer time span than only up to t=0.5. I thought that I could use a for loop for this purpose, but when I run my code I get the warning: Failure at t=5.086967e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.

Comment on my code:

- It is pretty much copy past of the example from the link above.

- For some reason, the example solving the DAE using ode15s in the link above do not use the result y0 returned by the decic function, instead the initial estimate y0est is used. If I try to simulate using the output from decic I get the error: Need a better guess y0 for consistent initial conditions. Which I find very strange since I thought decic would produce consistent initial conditions.

- I know that it's unnecessary to use ode15s two times, but for simplicity I perform the computation twice in order to get the plot automatically generated when no output arguments are given. I also know that it is unnecessary to have the code defining the equations inside the loop, but in another problem I'm working on I'm interested in implementing it this way. See code below:

clear all
close all
clc

m = 1.0;
r = 1.0;
g = 9.81;

y0est = [0.5*r; -0.8*r; 0; 0; 0; 0; 0];
yp0est = zeros(7,1);

results = zeros(7,6);
results(:,1) = y0est;

syms x(t) y(t) T(t);

for k=0:0.5:2

    vars = [x(t); y(t); T(t)];

    eqs= [m*diff(x(t), 2) == T(t)/r*x(t), ...
          m*diff(y(t), 2) == T(t)/r*y(t) - m*g, ...
          x(t)^2 + y(t)^2 == r^2];

    [eqs, vars, R] = reduceDifferentialOrder(eqs, vars);
    isLowIndexDAE(eqs,vars)
    [DAEs,DAEvars] = reduceDAEIndex(eqs,vars);
    [DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
    isLowIndexDAE(DAEs,DAEvars)

    [M,F] = massMatrixForm(DAEs,DAEvars);
    M = odeFunction(M, DAEvars);
    F = odeFunction(F, DAEvars);

    opt = odeset('Mass', M, 'InitialSlope', yp0est,...
                'RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));

    [y0, yp0] = decic(@(t,y,yp) M(t,y)*yp - F(t,y), 0, y0est, [], yp0est, [], opt);

    ode15s(F, [k, k+0.5], y0, opt);
    for j = 1:numel(DAEvars)
      S{j} = char(DAEvars(j));
    end
    legend(S, 'Location', 'Best')
    grid on
    hold on

    [Tresult,Yresult] = ode15s(F, [k, k+0.5], y0est, opt);
    results(:,(k+1)*2) = Yresult(end,:);

    y0est = Yresult(end,:);
end