function [x,t] = vanderPol(epsilon,tfinal)

%% The differential equation

f = @(t,x) [x(2), -x(1) - epsilon*(x(1).^2 - 1)*x(2)];
% systems of ODEs can be defined with function handles

%% Solving using regular Runge-Kutta (4th order)

x0 = [0.5, 0];       % initial conditions
[x,t] = RK4(f,[0 tfinal],x0);

plot(t,x(:,1),'-',t,x(:,2),'-.','LineWidth',2)
title('vanderPol oscillator','fontsize',18)
xlabel('time','fontsize',18)
legend('position','velocity')
figure
plot3(x(:,1),x(:,2),t,'LineWidth',2)
zlabel('time','fontsize',18)

xlabel('position','fontsize',18)
ylabel('velocity','fontsize',18)

figure
plot(x(:,1),x(:,2),'LineWidth',2)
xlabel('position','fontsize',18)
ylabel('velocity','fontsize',18)

end