function forced_vderpol

clear all;
close all;
clc;

%  ----------  Parameters  ----------
epsilon = 1.0; 
F       = 6.0;
w       = 2.0;
w1      = 2.0;

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

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


%  -----------  Numerical Solution  -------------
[t,v]=ode45(@vderpol_ODE,tspan,V0,[],epsilon,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 -5.0 5.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,[],F,w1);

rho   = vv(:,1); 
theta = vv(:,2);

figure(2);
subplot(2,1,1);
plot(t,rho,'b','LineWidth',2),
legend('\rho');
xlabel('Time','FontSize',16);
ylabel('\rho','FontSize',16);
axis([tmin tmax 0.0 4.0]);

subplot(2,1,2)
plot(rho,theta,'-');
xlabel('\rho','FontSize',16);
ylabel('\theta','FontSize',16);


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

u = rho.*cos(w*t-theta);

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



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

f=f';
    

%  ---------  Right-Hand-Side of Transformed ODEs  ---------- 
function g = transformed_ODE(t,vv,F,w1)    
rho   = vv(1);
theta = vv(2);
    
g(1) = rho*(4. - rho^2)/8. + (F/2.)*sin(theta);
g(2) = w1 + (F/(2.*rho))*cos(theta);

g=g';
