function cpendula (th1,th2,k,t,ti,d,m,l)
dt=t/ti; x1=(0:t/ti); x2=(0:t/ti); v1=(0:t/ti); v2=(0:t/ti);
x1(1)=th1; %set initial position to given thete coordinate x2(1)=th2;
thetaAccel1(theta1,theta2)=(k/(ml2))(d-sqrt(ll(cos(theta1)-cos(theta2))2))((l2(cos(theta1)-cos(theta2))sin(theta1))-lcos(theta1)(d+lsin(theta2)-lsin(theta1)))/sqrt(l2(cos(theta1)-cos(theta2)2)+(d+lsin(theta2)-lsin(theta1))2)-mgl*sin(theta1);
thetaAccel2(theta1,theta2)=(k/(ml2))(d-sqrt(ll(cos(theta1)-cos(theta2))2))((l2(cos(theta1)-cos(theta2))sin(theta2))-lcos(theta2)(d+lsin(theta2)-lsin(theta1)))/sqrt(l2(cos(theta1)-cos(theta2)2)+(d+lsin(theta2)-lsin(theta1))2)-mgl*sin(theta2);
%Need to have the inputs be initial theta, spring constant, initial %velocity.
for f= 2:ti;
v1(f)=v1(f-1)+thetaAccel1(x1(f-1),x2(f-1))*dt;
v2(f)=v2(f-1)+thetaAccel2(x1(f-1),x2(f-1))*dt;

x1(f)=x1(f-1)+v1(f-1)*dt+.5*thetaAccel1(x1(f-1),x2(f-1))*dt^2;
x2(f)=x2(f-1)+v2(f-1)*dt+.5*thetaAccel2(x1(f-1),x2(f-1))*dt^2;
end
plot(x1,dt);
end