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;
ZnVuY3Rpb24gZj12ZHBXSmFjb2JpYW4obXUsIEEsIHcsIHQsIFgpCgolUmVhcnJhbmdlIGlucHV0IGRhdGEKeD1YKDEpOyB5PVgoMik7ClE9W1goMyksIFgoNSk7CiAgIFgoNCksIFgoNildOwoKZj16ZXJvcyg2LDEpOwolRHVmZmluZydzIFZhbiBEZXIgUG9sIGVxdWF0aW9uCmYoMSk9eTsKZigyKT1tdSooMS14XjIpKnkteCtBKmNvcyh3KnQpOwoKJUxpbmVhcml6ZWQgc3lzdGVtIChKYWNvYmlhbikKSiA9IFsgICAgICAgICAgMCwgICAgICAgICAxOwogICAgLTIqbXUqeCp5LTEsIG11KigxLXheMildOwoKJVZhcmlhdGlvbmFsIGVxdWF0aW9uCmYoMzo2KT1KKlE7Cg==