function DAEpendulum()
%Simple pendulum using Diff. Alg. Eqs.
x0= 1; y0 = 0; vx0= 0; vy0= 0; % Initial conditions
z0= [ x0 y0 vx0 vy0]';         % Pack ICs into z0
tspan = linspace(0,10,101);    % Times to look at solution
p.m = 1; p.g = 1;              % Parameters

% ODESET sets preferences for ODE solvers like ODE45
options = odeset('abstol', 1e-6,'reltol',1e-6);

% The next line solves the odes that are in PENDRHS.
[tarray zarray] = ode45(@pendrhs,tspan, z0,options,p);

% Pull out x and y comps of position at different times
x= zarray(:,1);   y= zarray(:,2); t=tarray;

shg  % bring graphs to front
figure(1);  plot(t, x, t, y);  % Time history of x and y
figure (2); plot(x, y);        axis('equal')  % trajectory
figure (3)
%Animation in this loop
for i=1: length(t);      
    pause(.03);    %Pause for .03 s between movie frames
    plot(x(i), y(i), 'o','markersize',50)  % The pendulum is a dot.
     axis('equal')                  % make scale the same for x and y
     axis([ -1.2 1.2 -1.2 1.2]);    % Keep the frame with fixed axis  
end
end


%%%%%%%%%
function zdot = pendrhs(t,z,p)
m= p.m;  g = p.g;  % use simpler names for parameters
x=z(1); y= z(2); xdot = z(3);  ydot = z(4);  % unpack z into readable names
r = sqrt(x^2 + y^2);  v = sqrt(xdot^2 + ydot^2); %intermediate variables
A = [ m 0 x/r;     % Setting up the equation A*w = b
      0 m y/r;
      x y  0];
b = [0 -m*g -v^2]'; % The right hand side of the linear algebra equations

w = A\b;   %Backslash. The greatest thing in Matlab!
xddot = w(1);  yddot = w(2);  % pull out the accelerations and ignore the Tension
zdot = [xdot  ydot xddot yddot]';   % The rates of change of all variables in z
end