function [y1,m1,T,y0] = prob5(t,y,jaco)
%
% [alpha,m1,T,y0] = prob5
%
% or
%            f = prob5(t,y)
%
% or
%           f' = prob5(t,y,1)
%
% Problem 5 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
%
% Multi-term fractional version of:
%
% P.Panja. Dynamics of a fractional order predator-prey model with intra- 
% guild predation. Int. J. Math. Model. Simul. Appl. 39(4) (2019) 256–268. 
% https://doi.org/10.1080/02286203.2019.1611311         
%
% 
alfa1 = 0.99;
alfa2 = 0.8; 
y0    = [0.7 0.2 0.1];    
r     = [5 1 0.1];
a     = [0.01 1 35; 1 0.2 1; 0.1 1 0.3];
bet   = 0.01;
if nargin==0    % problem data
    y1    = [alfa1 alfa2];
    m1    = 1;
    T     = 100;
elseif nargin==1  
    y1    = [];
elseif nargin==2 % vector field
    if length(t)>1
       fatt  = 1./(1+bet*y(:,2));
       y1    = [r(1)*y(:,1)-a(1,1)*y(:,1).^2-a(1,2)*y(:,1).*y(:,2)-a(1,3)*y(:,1).*y(:,3) ...
             a(2,1)*y(:,1).*y(:,2)-a(2,2)*y(:,2).^2-a(2,3)*y(:,2).*y(:,3).*fatt-r(2)*y(:,2) ...
             a(3,1)*y(:,1).*y(:,3)+a(3,2)*y(:,2).*y(:,3).*fatt-a(3,3)*y(:,3).^2-r(3)*y(:,3)];
    else
       fatt  = 1/(1+bet*y(2));
       y1    = [r(1)*y(1)-a(1,1)*y(1)^2-a(1,2)*y(1)*y(2)-a(1,3)*y(1)*y(3); ...
             a(2,1)*y(1)*y(2)-a(2,2)*y(2)^2-a(2,3)*y(2)*y(3)*fatt-r(2)*y(2); ...
             a(3,1)*y(1)*y(3)+a(3,2)*y(2)*y(3)*fatt-a(3,3)*y(3)^2-r(3)*y(3)];
    end        
else             % Jacobian
    f1   = (1+bet*y(2));
    f2   = f1^2;
    y1   = [r(1)-2*a(1,1)*y(1)-a(1,2)*y(2)-a(1,3)*y(3)  -a(1,2)*y(1)  -a(1,3)*y(1); ...
            a(2,1)*y(2)  a(2,1)*y(1)-2*a(2,2)*y(2)-r(2)-a(2,3)*y(3)/f2  -a(2,3)*y(2)/f1; ...
            a(3,1)*y(3)  a(3,2)*y(3)/f2  a(3,1)*y(1)+a(3,2)*y(2)/f1-2*a(3,3)*y(3)-r(3)];
end
return