function [y1,T,y0] = ex4(t,y,jaco)
%
% [alpha,T,y0] = ex4
%
% or
%            f = ex4(t,y)
%
% or
%           f' = ex4(t,y,1)
%
% or
%           yT = ex4(5)
%
% Example 4 in:
%
%  L.Brugnano, G.Gurioli, F.Iavernaro. Numerical solution of FDE-IVPs by 
%  using fractional HBVMs: the fhbvm code. Numerical Algorithms (2024). 
%  https://doi.org/10.1007/s11075-024-01884-y
%
alfa = 0.7;
A    = 1;
B    = 3;
if nargin==0    % problem data
    y1   = alfa;
    T    = 5;
    y0   = [1.2 2.8];
elseif nargin==1 % y(5), numerically evaluated with FHBVM(24,22), n=1000,
    y1   = [8.904632063462272e-01  3.326603532694057e+00];     %  h0=1e-14.
elseif nargin==2 % vector field
    if length(t)>1
       y1 = [A-(B+1)*y(:,1)+y(:,1).^2.*y(:,2) B*y(:,1)-y(:,1).^2.*y(:,2)];
    else
       y1 = [A-(B+1)*y(1)+y(1)^2*y(2); B*y(1)-y(1)^2*y(2)];
    end        
else             % Jacobian
    y1   = [-(B+1)+2*y(1)*y(2) y(1)^2; B-2*y(1)*y(2) -y(1)^2];
end
return