function [t,y]=Double_Pendulum(y10,y20,y30,y40);

%%Double_Pendulum(phi1(0),phi1_dot(0),phi2(0),phi2_dot(0)--This function will take in the initial conditions 
%%for the Nonlinear double pendulum system and plot the positions of the
%%masses

%%Representing Cos(theta2-theta1) as c1
%              Sin(theta2-theta1) as s1


y1='[y(2);';
y2='s1*(y(4)^2+(y(2))^2*c1)+9.8*(sin(y(3))*c1-2*sin(y(1)))/(2-(c1)^2);';
y3='y(4);';
y4='s1*(c1*y(4)^2+2*y(2)^2)-2*9.8*(sin(y(1)*c1-sin(y(3)))/((c1)^2-2))]';

F=strcat(y1,y2,y3,y4);

c1='(cos(y(3)-y(1)))';
s1='(sin(y(3)-y(1)))';


F_a=regexprep(F,'c1',c1);
F_b=regexprep(F_a,'s1',s1);

%if length(findobj)==1
 %   i=0;
 %else
 %   delete(findobj(1));
  %  delete(findobj(2));
  %end

markers=['--' '.' '*'];
func=inline(F_b,'t','y');

[t,y]=ode45(func,[0.01,10],[y10,y20,y30,y40]);

clf reset;

subplot(2,2,1);
plot(y(:,1),y(:,3));
xlabel('\theta_1');
ylabel('\theta_2');

subplot(2,2,2);
plot(t,y(:,1));
xlabel('t');ylabel('\theta_1');

subplot(2,2,3);
plot(t,y(:,3));
xlabel('t');ylabel('\theta_2');