%%%Dual Harmonic oscillators
%2 frequencies
clear all
close all
w1 = 1;
w2 = (1+sqrt(5))/2;

w1s = w1^2;
w2s = w2^2;

ff = @(t,X)  [X(2); -w1s.*X(1); X(4); -w2s.*X(3)];


tf = 1000;
T = [0 tf];

%IC
X0 = [1,0,0,1];

%%Solve the ode
options = odeset('RelTol',1e-16,'AbsTol',1e-16,'MaxStep',1e-7)
[tt,XX] = ode45(ff,T,X0);

%plot positions
figure(1)
plot(XX(:,1),XX(:,3))


figure(2)
subplot(2,1,1)
plot(XX(:,1),XX(:,2))
subplot(2,1,2)
plot(XX(:,3),XX(:,4))
xlim([-1.5 1.5])
ylim([-1.5 1.5])


%Exact solution
TT = [0:0.002:tf];
x1 = X0(1)*cos(w1*TT);
x2 = X0(4)*sin(w2*TT);

figure(3)
plot(x1,x2)