function [y1,T,y0] = vdp(t,y,jaco)
%
% Fractional Van der Pol oscillator
%
% [alpha,T,y0] = vdp
%
% or
%           f = vdp(t,y)     -- vector field (also in vector mode)
%
% or
%           f' = vdp(t,y,1)  -- Jacobian
%
% or
%            y = vdp(t)      -- reference solution at t=30
%
% Example 3 in:
%
% L.Brugnano, G.Gurioli, F.Iavernaro, M.Vikerpuur. 
% Analysis  and  implementation of  collocation methods for
% fractional differential equations.
% J. Sci. Comput. 104 (2025) 92.
% https://doi.org/10.1007/s10915-025-03006-9
%
alfa = 0.9;
epsi = 10;
y00  = [0 -2];
if nargin==0    % problem data
    y1   = alfa;
    T    = 30;
    y0   = y00;
elseif nargin==1 % reference solution at t=30
    %computed by: [t,y,et]=fhbvm2('vdp',y0,T,10000,5,100);
    y1  = [-9.731412336030406e-01     4.817675882545460e-01];
elseif nargin==2 % vector field
    if length(t)>1
      y1 = [y(:,2) -y(:,1)-epsi*( y(:,1).^2-1 ).*y(:,2) ];
    else
      y1 = [y(1); -y(1)-epsi*(y(1)^2-1)*y(2)];  % a column vector
    end  
else             % Jacobian
    y1   = [0 1;-1-2*epsi*y(1)*y(2) -epsi*(y(1)^2-1)];
end
return
end