function [y1,T,y0] = exe3(t,y,jaco)
%
%  [alpha,T,y0] = exe3
% 
%  or 
%             f = exe3(t,y)
% 
%  or 
%            f' = exe3(t,y,1)
% 
%  or
%             y = exe3(t)
%
%  Example 3 in   https://arxiv.org/abs/2407.11460
%
%  Adapted from Examples 2 and 3 in
%
%  T.B.Nguyen, B.Jang. A high-order predictor-corrector method for 
%    solving nonlinear differential equations of fractional order. 
%    Fractional Calculus and Applied Analysis 20(2) (2017) 447–476. 
%    https://doi.org/10.1515/fca-2017-0023
%
alfa = 1.25;
if nargin==0    % problem data
    y1   = alfa;
    T    = 1;
    y0   = [0 0;0 0];
elseif nargin==1 % solution
    y1   = [t.^(3+alfa) t.^(4+alfa)];
elseif nargin==2 % vector field
%    y1   = [gamma(4+alfa)/6*t.^3-t.^(8+2*alfa)+y(:,2).^2   ...
%            gamma(5+alfa)/4*t.^4+t.^(3+alfa)-y(:,1)];
   if length(t)>1    

      y1 = [gamma(4+alfa)/6*t.^3-t.^(8+2*alfa)+y(:,2).^2   ...
            gamma(5+alfa)/24*t.^4+t.^(3+alfa)-y(:,1)];
    
   else
   
      y1 = [gamma(4+alfa)/6*t^3-t^(8+2*alfa)+y(2)^2;   
            gamma(5+alfa)/24*t^4+t^(3+alfa)-y(1)];    
       
   end    
else             % Jacobian
    y1   = [0 2*y(2); -1 0]; 
end
return