function run_vdpol_linear
%% Linearised Van der Pol's oscillator
%
% $$
% x_1' = x_2
% x_2' = -x_1 + 2 a x_2
% $$

a = -0.1;
figure;

%% Initial condition $x_0=(1,1)^T$

x0 = [1;1];
% Plot initial point
plot(x0(1),x0(2),'r*');
hold on;
% Solve forward
[~,x]=ode23s(@(t,x)vdpol_linear(t,x,a), [0,20], x0);
plot(x(:,1),x(:,2),'r');
% Solve backward
[~,x]=ode23s(@(t,x)vdpol_linear(t,x,a), [0,-1.5], x0);
plot(x(:,1),x(:,2),':r');

%% Initial condition $x_0=(1.5,1)^T$

x0 = [1.5;1];
% Plot initial point
plot(x0(1),x0(2),'b*');
hold on;
% Solve forward
[~,x]=ode23s(@(t,x)vdpol_linear(t,x,a), [0,200], x0);
plot(x(:,1),x(:,2),'b');
% Solve backward
[~,x]=ode23s(@(t,x)vdpol_linear(t,x,a), [0,-1.5], x0);
plot(x(:,1),x(:,2),':b');

%% Initial condition $x_0=(0.5,1)^T$

x0 = [0.5;1];
% Plot initial point
plot(x0(1),x0(2),'m*');
hold on;
% Solve forward
[~,x]=ode23s(@(t,x)vdpol_linear(t,x,a), [0,200], x0);
plot(x(:,1),x(:,2),'m');
% Solve backward
[~,x]=ode23s(@(t,x)vdpol_linear(t,x,a), [0,-1.5], x0);
plot(x(:,1),x(:,2),':m');

xlabel('x_1');
ylabel('x_2');
xlim([-2 2]);
ylim([-2 4]);