function [y1,m1,T,y0] = prob2(t,y,jaco)
%
% [alpha,m1,T,y0] = prob2
%
% or
%            f = prob2(t,y)
%
% or
%           f' = prob2(t,y,1)
%
% or
%            y = prob2(t)
%
% Problem 2 in:
%
%  L.Brugnano, G.Gurioli, F.Iavernaro, M.Vikerpuur. 
%  A Multi-Order Extension of Fractional HBVMs (FHBVMs)
%  J. Sci. Comput. 107 (2026) 94. https://doi.org/10.1007/s10915-026-03315-7
%
%
% Adapted from Problem 7 in:
%
%  L.Brugnano, G.Gurioli, F.Iavernaro, M.Vikerpuur. 
%  FDE-testset: comparing Matlab© codes for solving fractional
%  differential equations of Caputo type
%  Fractal Fract. 2025, 9(5), 312; https://doi.org/10.3390/fractalfract9050312
%
% 
alfa1 = 0.2;
alfa2 = 0.4;
beta0 = 0.1;

sol = @(t,alfa,bet) (1-t.^2).^2 +4*t.^alfa ...
                  + (2-3*t.^(0.2)).*t.^(alfa+bet); 

g   = @(t,alfa,bet) ( 24/gamma(5-alfa) )*t.^(4-alfa) ...
                  - ( 4/gamma(3-alfa) )*t.^(2-alfa) ...
                - ( 3*gamma(1.2+alfa+bet)/gamma(1.2+bet) )*t.^(0.2+bet) ...
                + ( 2*gamma(1+alfa+bet)/gamma(1+bet) )*t.^bet ...
                + 4*gamma(1+alfa);
if nargin==0    % problem data
    y1    = [alfa1 alfa2];
    m1    = 1;
    T     = 2;
    y0    = [sol(0,alfa1,beta0) sol(0,alfa2,beta0)];
elseif nargin==1 % solution
    y1    = [sol(t,alfa1,beta0) sol(t,alfa2,beta0)];
elseif nargin==2 % vector field
    if length(t)>1
       s1 = sol(t,alfa1,beta0);
       s2 = sol(t,alfa2,beta0);
       y1 = [s2.^2-y(:,2).^2+g(t,alfa1,beta0) ...
             -s1.^2+y(:,1).^2+g(t,alfa2,beta0)];
    else
       s1 = sol(t,alfa1,beta0)';
       s2 = sol(t,alfa2,beta0)';
       y1 = [s2^2-y(2)^2+g(t,alfa1,beta0)'; ...
             -s1^2+y(1)^2+g(t,alfa2,beta0)'];
    end    
else             % Jacobian
    y1    = [0 -2*y(2); 2*y(1) 0];  
end
return