%Perturbation_VanderPol_fun.m
function dy=Perturbation_VanderPol_fun(t,y,epsilon)

y1=y(1);
y2=y(2);

dy1=y2;
dy2=-y1+epsilon*(1-y1^2)*y2;

dy=[dy1; dy2];
