function vderpol

clear all;
close all;
clc;

%  ----------  Parameters  ----------
mu   = 1.0;
tmin = 0.0;  tmax = 50;
xmin = -2.5; xmax = 2.5; 
h    = 0.01;

trans       = 1.0;
tspan_trans = [tmin:h:trans];
tspan       = [tmin:h:tmax];

%  ----------  Intial Conditions ----------
V0   = [0.26, 0.1]; 
[t_trans,v_trans] = ode45(@rhs_ode,tspan_trans,V0,[],mu);


%  -----------  Integrate ODEs  -----------
V0    = v_trans(end,:)';
[t,v]=ode45(@rhs_ode,tspan,V0,[],mu);


%  -----------  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]);

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


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

f=f';
    