function f=vdpWJacobian(mu, A, w, t, X)

%Rearrange input data
x=X(1); y=X(2);
Q=[X(3), X(5);
   X(4), X(6)];

f=zeros(6,1);
%Duffing's Van Der Pol equation
f(1)=y;
f(2)=mu*(1-x^2)*y-x+A*cos(w*t);

%Linearized system (Jacobian)
J = [          0,         1;
    -2*mu*x*y-1, mu*(1-x^2)];

%Variational equation
f(3:6)=J*Q;