function [y1,T,y0] = ex3(t,y,jaco)
%
% [alpha,T,y0] = ex3
%
% or
%            f = ex3(t,y)
%
% or
%           f' = ex3(t,y,1)
%
% or
%            y = ex3(t)
%
% Example 3 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 = 1/3;
if nargin==0    % problem data
    y1   = alfa;
    T    = 1;
    y0   = [1 0];
elseif nargin==1 % solution
    y1   = [t.^(2/3)+1 t.^(4/3)];
elseif nargin==2 % vector field
%    y11  = t.*y(:,1).^3/10 - t.*(t.^(2/3)+1).^3/10 + ( gamma(5/3)/gamma(4/3) )*t.^(1/3);
%    y12  = y(:,2).^3/3 - t.^4/3 + gamma(7/3)*t;
    if length(t)>1
       y11  = t.*y(:,1).^3/10 - t.*(sqrt(abs(y(:,2)))+1).^3/10 ...
              + ( gamma(5/3)/gamma(4/3) )*t.^(1/3);
       y12  = y(:,2).^3/3 - (y(:,1)-1).^6/3 + gamma(7/3)*t;
       y1   = [y11 y12];
    else
       y11  = t*y(1)^3/10 - t*(sqrt(abs(y(2)))+1)^3/10 ...
              + ( gamma(5/3)/gamma(4/3) )*t^(1/3);
       y12  = y(2)^3/3 - (y(1)-1)^6/3 + gamma(7/3)*t;
       y1   = [y11; y12];
    end    
else  % Jacobian
    y1  = [(3/10)*t*y(1)^(3/10-1) (3/10)*t*sqrt(abs(y(2))+1); 
           -2*(y(1)-1)    y(2)^2];
end
return