function [ dy ] = harmonic( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1))/m;
end


function [ dy ] = anharmonicA( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1)+l*y(1)^2)/m;
end



function [ dy ] = anharmonicB( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1)+l*y(1)^3)/m;
end





m=1;
k=1;
l=0.2;
c=0;

x0=1;
x01=2.4;
x02=3;
p0=0;



dy=@(t,y) anharmonicA(t,y,m,k,l,c);
dyB=@(t,y) anharmonicB(t,y,m,k,l,c);
dyh=@(t,y) harmonic(t,y,m,k,l,c);
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);

PSRange=[-10 10];

subplot(3,2,1);
[T,Y] = ode45(dy,[0 40],[x0, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x0, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x0, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
ylim([-1.5 1.5]*x0);
xlabel('time t');
ylabel('position x(t)');
title(['solution for x_0=' num2str(x0) ', p_0=0']);

subplot(3,2,2);
for x0=-5:0.5:10
    [Th,Yh] = ode45(dyh,[0 100],[x0, p0],options);
    plot(Yh(:,1),Yh(:,2), 'b-');
    if x0==0.5
        hold on
    end
end
hold off
xlim(PSRange);
ylim(PSRange);
title(['phase space plot, harmonic oscillator']);


subplot(3,2,3);
[T,Y] = ode45(dy,[0 40],[x01, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x01, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x01, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
xlabel('time t');
ylabel('position x(t)');
title(['solution for x_0=' num2str(x01) ', p_0=0']);
ylim([-2 2]*x01);
subplot(3,2,4);
for x0=-5:0.5:10
    [T,Y] = ode45(dy,[0 100],[x0, p0],options);
    plot(Y(:,1),Y(:,2), 'r-');
    if x0==0.5
        hold on
    end
end
hold off
xlim(PSRange);
ylim(PSRange);
title(['phase space plot, asymmetric force law']);


subplot(3,2,5);
[T,Y] = ode45(dy,[0 40],[x02, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x02, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x02, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
xlabel('time t');
ylabel('position x(t)');
title(['solution for x_0=' num2str(x02) ', p_0=0']);
ylim([-3 3]*abs(x02));
subplot(3,2,6);
for x0=-5:0.5:10
    [TB,YB] = ode45(dyB,[0 100],[x0, p0],options);
    plot(YB(:,1),YB(:,2), 'g-');
    if x0==0.5
        hold on
    end
end
hold off
xlim(PSRange);
ylim(PSRange);
xlabel('position x');
ylabel('impulse p');
title(['phase space plot, symmetric force law']);