function [y1,T,y0] = exe1(t,y,jaco)
%
%  [alpha,T,y0] = exe1
% 
%  or
%             f = exe1(t,y)
% 
%  or
%            f' = exe1(t,y,1)
% 
%  or
%             y = exe1(t)
%
%  Example 1 in https://arxiv.org/abs/2407.11460
%
%  Taken from Example C.1, pag. 207 in:
%
%  K.Diethelm. The analysis of fractional differential equations. 
%  An application-oriented exposition using differential operators of 
%  Caputo type. Lecture Notes in Math. 2004. Springer-Verlag, Berlin, 2010.
%
alfa = 1.3;
if nargin==0    % problem data
    y1   = alfa;
    T    = 1;
    y0   = [0;0];
elseif nargin==1 % solution
    y1   = t.^8-3*t.^(4+alfa/2)+(9/4)*t.^alfa;
elseif nargin==2 % vector field
    y1   = -abs(y).^(3/2) +( 40320/gamma(9-alfa) )*t.^(8-alfa) -...
           3*( gamma(5+alfa/2)/gamma(5-alfa/2) )*t.^(4-alfa/2) +...
           ( (3/2)*t.^(alfa/2) - t.^4 ).^3 + (9/4)*gamma(alfa+1);
else             % Jacobian
    y1   = -1.5*sqrt(abs(y));
end
return