function forced_duffing

clear all;
close all;
clc;

%  ----------  Parameters  ----------
epsilon = 0.1; 
c       = 0.0;
alpha   = 10.0;
F       = 4.0;
w       = 8.0;
w1      = 8.0;

%  ----------  Intial Conditions ----------
V0  = [0.26, 0.1];
VV0 = [2.0, 1.0]; 

%  -----------  Integrate ODEs  -----------
tmin = 0.0;  
tmax = 100;
xmin = -4.0; 
xmax = 4.0;
h    = 0.05;
tspan = [tmin:h:tmax];


%  -----------  Numerical Solution  -------------
[t,v]=ode45(@duffing_ODE,tspan,V0,[],epsilon,c,alpha,F,w);

x = v(:,1); 
y = v(:,2);

figure(1);
subplot(2,1,1);
plot(t,x,'b',t,y,'r--','LineWidth',2),
legend('x','y');
xlabel('Time','FontSize',16);
ylabel('x,y','FontSize',16);
title('Numerical Solution');
axis([tmin tmax -1.0 1.0]);

subplot(2,1,2)
plot(x,y);
xlabel('x','FontSize',16);
ylabel('y','FontSize',16);
 

%  -----------  Averaging ODE  -------------
[t,vv]=ode45(@transformed_ODE,tspan,VV0,[],c,w1,alpha,F);

A = vv(:,1); 
B = vv(:,2);

figure(2);
subplot(2,1,1);
plot(t,A,'b',t,B,'r--','LineWidth',2),
legend('A','B');
xlabel('Time','FontSize',16);
ylabel('A,B','FontSize',16);
axis([tmin tmax -4.0 4.0]);

subplot(2,1,2)
plot(A,B,'-');
xlabel('A','FontSize',16);
ylabel('B','FontSize',16);


%  -----------  Multiple Scale Solution  -------------
figure(3);

u = A.*cos(w*t) + B.*sin(w*t);

plot(t,u,'b','LineWidth',2),
legend('u_{MS}');
xlabel('Time','FontSize',16);
ylabel('u_{(t;\epsilon)}','FontSize',16);
title('Multiple Scale Solution');
axis([tmin tmax xmin xmax]);



%  ---------  Right-Hand-Side of Forced Duffing ODEs  ---------- 
function f = duffing_ODE(t,v,epsilon,c,alpha,F,w)    
x    = v(1);
y    = v(2);
    
f(1) = y;
f(2) = -x - epsilon*c*y - epsilon*alpha*x^3 + epsilon*F*cos(w*t);

f=f';
    

%  ---------  Right-Hand-Side of Transformed ODEs  ---------- 
function g = transformed_ODE(t,vv,c,w1,alpha,F)    
A   = vv(1);
B   = vv(2);
Rs  = A^2 + B^2;
    
g(1) = -(c/2.)*A - w1*B    + (3./8.)*alpha*B*Rs;
g(2) = w1*A      -(c/2.)*B - (3./8.)*alpha*A*Rs + F;

g=g';
