%Lecture35.m

clear all

m=input('Enter m: ');
k=input('Enter k: ');
I1=m;
I2=2/3*m;
C=1/k;

p=[I1 I2 C];

x0(1)=input('\nInitial Condition for q1: ');
x0(2)=input('Initial Condition for p1: ');
x0(3)=input('Initial Condition for q2: ');
x0(4)=input('Initial Condition for p2: ');

tf=20;
[t,x]=ode45('Lecture35_Equations',[0 tf],x0,[],p);


figure(1)
plot(t,x(:,1),'r',t,x(:,3),'b:')
legend('q_{1}','q_{2}')

figure(2)
plot(t,x(:,1),'r',t,x(:,1)-x(:,3),'b:')
legend('x_{1}','x_{2}')


h=0.1;
N=tf/h;
figure(3)
plot([-2 -2],[-1 1],'k:',[2 2],[-1 1],'k:')
hold on
plot(-2+x0(1),0,'ro',2+x0(1)-x0(3),0,'bo')
hold off
axis([-5 5 -1 1])
pause

for i=1:N+1
   time=(i-1)*h;
   
   q1=spline(t,x(:,1),time);
   q2=spline(t,x(:,3),time);
      
	plot([-2 -2],[-1 1],'k:',[2 2],[-1 1],'k:')
   hold on
   plot(-2+q1,0,'ro',2+q1-q2,0,'bo')
   hold off
   
   axis([-5 5 -1 1])
   pause(0.1)
end

