function vderpol

clear all;
close all;
clc;

%  ----------  Parameters  ----------
epsilon = 1.0; 
a       = 0.75;

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

%  -----------  Integrate ODEs  -----------
tmin = 0.0;  
tmax = 30;
xmin = -2.5; 
xmax = 2.5;
h    = 0.1;
tspan = [tmin:h:tmax];

[t,v]=ode45(@vderpoleq,tspan,V0,[],epsilon);

%  -----------  Time Series Solution  -------------
x = v(:,1); 
y = v(:,2);

figure(1);
subplot(2,1,1);
plot(t,x,'r','LineWidth',2),  ylabel('x','FontSize',16);
axis([tmin tmax xmin xmax]);
  
subplot(2,1,2);
plot(t,y,'b','LineWidth',2), ylabel('y','FontSize',16);
axis([tmin tmax xmin xmax]);

%  ----------  Assymptotic Solution  ----------
figure(2);
u_appx = cos(t)- epsilon*((a^3-4.0*a)/8.)*t.*cos(t) - ...
         (epsilon/32.)*a^3*sin(3*t);
plot(t,x,t,u_appx,'LineWidth',2);
legend('x','u_{appx}');

%  -----------  Phase Plane  -------------
figure(3)
plot(x,y);
xlabel('x','FontSize',16);
ylabel('y','FontSize',16);


%  ---------  Right-Hand-Side of ODEs  ---------- 
function f = vderpoleq(t,v,epsilon)    
x    = v(1);
y    = v(2);
    
f(1) = y;
f(2) = -x + epsilon*(1.0-x^2)*y;

f=f';
    