function  twomassesconstrained()
%Two masses connected by a Rod
%Example from lecture on 9/10/2012
%NOT ADEQUATELY COMMENTED
%parameters
p.F0 = 1;  p.m = [1 1]'; p.L0=3; p.F0=2;
z0=[0 p.L0 0 0]';
tspan = linspace(0,100,1001);

%NOTE NEW FEATURE: ODESET
options = odeset('reltol',1e-13,'abstol',1e-13);
[t zarray] = ode23(@rhs, tspan,z0, options,p);

x1=zarray(:,1);
x2=zarray(:,2);

subplot(2,1,1)
plot(t,x1,'r', t,x2,'b')
ylabel('Position')

subplot(2,1,2)
plot(t,x2-x1-p.L0)
axis([ -inf inf -1e-13 1e-13])
ylabel ('constraint error')

end



function zdot = rhs(t,z,p)
if t >= 0
    F = p.F0;
end
x= z(1:2);
v= z(3:4);

M=[p.m(1) 0; 0 p.m(2)];
J= [-1 1]';
N = [M  J; 
     J' 0];     % Implicitly holds left side of ODEs and Const Eqs.
b = [ 0 F 0]';  % right hand side of linear equations

%THIS IS THE KEY LINE, Solving for accelerations
w =N\b; %  the list of accelerations and the constraint Tension

xdot = v;        % both masses
vdot = w(1:2);   % w(3) is T, but we have no need for it

zdot = [xdot; vdot];
    
end