function rhs = Fvdp(t,y)

%
%  save file as Fvdp.m
%     - right-side forcing for ode45 solver
%

global ODE ep mu

rhs = zeros(size(y));

switch ODE
    case{1}
        rhs(1) =  y(2);
        rhs(2) = -mu*(y(1)^2 - 1)*y(2) - y(1);
    case{2}
        rhs(1) =  mu*(y(2) - (y(1)^2 - 3)*y(1)/3);
        rhs(2) = -y(1)/mu;
end