%% Solve Differential Algebraic Equations (DAEs)
%% 1

syms x(t) y(t) T(t) m r g

eqn1 = m*diff(x(t), 2) == T(t)/r*x(t);
eqn2 = m*diff(y(t), 2) == T(t)/r*y(t) - m*g;
eqn3 = x(t)^2 + y(t)^2 == r^2;
eqns = [eqn1 eqn2 eqn3];

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

%
%% 2

M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);

isLowIndexDAE(eqns,vars);

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);

isLowIndexDAE(DAEs,DAEvars);

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars);

f = daeFunction(DAEs, DAEvars, g, m, r);

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

F = @(t, Y, YP) f(t, Y, YP, g, m, r);

DAEvars;

y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0; 0; 0];
yp0est = zeros(7,1);

opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));

[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

%
%% 3

[tSol,ySol] = ode15i(F, [0, 0.5], y0, yp0, opt);
plot(tSol,ySol(:,1:origVars),'-o')

for k = 1:origVars
  S{k} = char(DAEvars(k));
end

legend(S, 'Location', 'Best')
grid on

%
%% 4

r = 2;
F = @(t, Y, YP) f(t, Y, YP, g, m, r);

y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0; 0; 0];
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

%
%% 5

[tSol,y] = ode15i(F, [0, 0.5], y0, yp0, opt);
plot(tSol,y(:,1:origVars),'-o')

for k = 1:origVars
  S{k} = char(DAEvars(k));
end
legend(S, 'Location', 'Best')
grid on

% Copyright 2012 The MathWorks, Inc.