global e;

e = 2;
x0 = 3*(rand(2,1)-1/2);
N = 300;
Tfin = 20;
time = linspace(0,Tfin,N);
[t,phi] = ode45(@vdpfun,time,x0);
r = 4;
clf

for k = 1:N
    plot(phi(:,1),phi(:,2),'b-',phi(k,1),phi(k,2),'ro')
    axis square
    axis([-r r -r r])
    pause(3/N)
end

pause(3)

plot(time,phi(:,1),'b-',time,phi(:,2),'r-')
