function [y1,T,y0] = prob1(t,y,jaco)
%
% [alpha,T,y0] = prob1
%
% or
%            f = prob1(t,y)
%
% or
%           f' = prob1(t,y,1)
%
% or
%            y = prob1(t)
%
% Problem 1 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 4 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
%
% 
alfa = 0.5;   % must be 0.5 %      
cc   = sqrt(2)/2;
ii   = sqrt(-1);
ll   = 10;
A0   = [1 1;-1 1];
Q    = eye(5); Q(2:3,2:3) = A0; Q(4:5,4:5) = A0;
Q1   = eye(5); Q1(2:3,2:3) = A0'/2; Q1(4:5,4:5) = A0'/2;
QQ   = eye(5); QQ(1:2,1:2) = A0; QQ(3:4,3:4) = A0;
QQ1  = eye(5); QQ1(1:2,1:2) = A0'/2; QQ1(3:4,3:4) = A0'/2;
V    = [[-ii ii; 1 1] zeros(2,3); ...
        zeros(2,2) [-ii ii; 1 1] zeros(2,1); [0 0 0 0 1];];
A    = QQ*Q*[ll*A0 zeros(2,3); zeros(2,2) A0/2 zeros(2,1); [0 0 0 0 -1];]*Q1*QQ1;
lam  = [ll*[1+ii 1-ii] [1+ii 1-ii]/2 -1];   % A = QQ*Q*V*diag(lam)/V*Q1*QQ1
y00  = [1 2 3 4 5];
if nargin==0        % problem data
    y1   = alfa;
    T    = 20;
    y0   = y00;
elseif nargin==1    % solution
    talf = t(:).^alfa;
    y1   = V\(Q1*(QQ1*y00(:)));
    yt1  = ml(lam(1)*talf,alfa)*y1(1);
    yt2  = ml(lam(2)*talf,alfa)*y1(2);
    yt3  = ml(lam(3)*talf,alfa)*y1(3);
    yt4  = ml(lam(4)*talf,alfa)*y1(4);
    yt5  = ml(lam(5)*talf,alfa)*y1(5);
    y1   = real([yt1 yt2 yt3 yt4 yt5]*V.'*Q.'*QQ.');
elseif nargin==2 % vector field
    if length(t)>1
      y1 = y*A';
    else
      y1 = A*y(:);  % a column vector
    end  
else                % Jacobian
    y1   = A;
end
return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% auxiliary function

function E = ml(z,alpha,beta,gama)
%
% Evaluation of the Mittag-Leffler (ML) function with 1, 2 or 3 parameters
% by means of the OPC algorithm [1]. The routine evaluates an approximation
% Et of the ML function E such that |E-Et|/(1+|E|) approx 1.0e-15   
%     
%
% E = ML(z,alpha) evaluates the ML function with one parameter alpha for
% the corresponding elements of z; alpha must be a real and positive
% scalar. The one parameter ML function is defined as
%
% E = sum_{k=0}^{infty} z^k/Gamma(alpha*k+1)
%
% with Gamma the Euler's gamma function.
%
%
% E = ML(z,alpha,beta) evaluates the ML function with two parameters alpha
% and beta for the corresponding elements of z; alpha must be a real and
% positive scalar and beta a real scalar. The two parameters ML function is
% defined as
%
% E = sum_{k=0}^{infty} z^k/Gamma(alpha*k+beta)
%
%
% E = ML(z,alpha,beta,gama) evaluates the ML function with three parameters
% alpha, beta and gama for the corresponding elements of z; alpha must be a
% real scalar such that 0<alpha<1, beta any real scalar and gama a real and
% positive scalar; the arguments z must satisfy |Arg(z)| > alpha*pi. The
% three parameters ML function is defined as 
%
% E = sum_{k=0}^{infty} Gamma(gama+k)*z^k/Gamma(gama)/k!/Gamma(alpha*k+beta)
%
%
% NOTE: 
% This routine implements the optimal parabolic contour (OPC) algorithm
% described in [1] and based on the inversion of the Laplace transform on a
% parabolic contour suitably choosen in one of the regions of analyticity
% of the Laplace transform.
%
%
% REFERENCES
%
%   [1] R. Garrappa, Numerical evaluation of two and three parameter
%   Mittag-Leffler functions, SIAM Journal of Numerical Analysis, 2015,
%   53(3), 1350-1369
%
%
%   Please, report any problem or comment to : 
%        roberto dot garrappa at uniba dot it
%
%   Copyright (c) 2015, Roberto Garrappa, University of Bari, Italy
%   roberto dot garrappa at uniba dot it
%   Homepage: http://www.dm.uniba.it/Members/garrappa
%   Revision: 1.4 - Date: October 8 2015


% Check inputs
if nargin < 4
    gama = 1 ;
    if nargin < 3
        beta = 1 ;
        if nargin < 2
            error('MATLAB:ml:NumberParameters', ...
                'The parameter ALPHA must be specified.');
        end
    end
end

% Check whether the parameters ALPHA, BETA and GAMMA are scalars
if length(alpha) > 1 || length(beta) > 1 || length(gama) > 1
    alpha = alpha(1) ; beta = beta(1) ; gama = gama(1) ;
    warning('MATLAB:ml:ScalarParameters', ...
        ['ALPHA, BETA and GAMA must be scalar parameters. ', ...
        'Only the first values ALPHA=%f BETA=%f and GAMA=%f will be used. '], ...
        alpha, beta, gama) ;
end

% Check whether the parameters meet the contraints
if real(alpha) <= 0 || real(gama) <= 0 || ~isreal(alpha) || ...
        ~isreal(beta) || ~isreal(gama)
    error('MATLAB:ml:ParametersOutOfRange', ...
        ['Error in the parameters of the Mittag-Leffler function. ', ...
        'Parameters ALPHA and GAMA must be real and positive. ', ...
        'The parameter BETA must be real.']) ;
end

% Check parameters and arguments for the three parameter case
if abs(gama-1) > eps
    if alpha > 1
        error('MATLAB:ml:ALPHAOutOfRange',...
            ['With the three parameters Mittag-Leffler function ', ...
            'the parameter ALPHA must satisfy 0 < ALPHA < 1']) ;
    end
    if min(abs(angle(z(abs(z)>eps)))) <= alpha*pi 
        error('MATLAB:ml:ThreeParametersArgument', ...
            ['With the three parameters Mittag-Leffler function ', ...
            'this code works only when |Arg(z)|>alpha*pi.']);
    end
end

% Target precision 
log_epsilon = log(10^(-15)) ; 

% Inversion of the LT for each element of z
E = zeros(size(z)) ;  
for k = 1 : length(z)
    if abs(z(k)) < 1.0e-15
        E(k) = 1/gamma(beta) ; 
    else
        E(k) = LTInversion(1,z(k),alpha,beta,gama,log_epsilon) ;
    end
end


end


% =========================================================================
% Evaluation of the ML function by Laplace transform inversion
% =========================================================================
function E = LTInversion(t,lambda,alpha,beta,gama,log_epsilon)

% Evaluation of the relevant poles
theta = angle(lambda) ;
kmin = ceil(-alpha/2 - theta/2/pi) ;
kmax = floor(alpha/2 - theta/2/pi) ;
k_vett = kmin : kmax ;
s_star = abs(lambda)^(1/alpha) * exp(1i*(theta+2*k_vett*pi)/alpha) ;

% Evaluation of phi(s_star) for each pole
phi_s_star = (real(s_star)+abs(s_star))/2 ;

% Sorting of the poles according to the value of phi(s_star)
[phi_s_star , index_s_star ] = sort(phi_s_star) ;
s_star = s_star(index_s_star) ;

% Deleting possible poles with phi_s_star=0
index_save = phi_s_star > 1.0e-15 ;
s_star = s_star(index_save) ;
phi_s_star = phi_s_star(index_save) ;

% Inserting the origin in the set of the singularities
s_star = [0, s_star] ;
phi_s_star = [0, phi_s_star] ;
J1 = length(s_star) ; J = J1 - 1 ;

% Strength of the singularities
p = [ max(0,-2*(alpha*gama-beta+1)) , ones(1,J)*gama ]  ;
q = [ ones(1,J)*gama , +Inf] ;
phi_s_star = [phi_s_star, +Inf] ;

% Looking for the admissible regions with respect to round-off errors
admissible_regions = find( ...
    (phi_s_star(1:end-1) < (log_epsilon - log(eps))/t) & ...
    (phi_s_star(1:end-1) < phi_s_star(2:end))) ;

% Initializing vectors for optimal parameters
JJ1 = admissible_regions(end) ;
mu_vett = ones(1,JJ1)*Inf ;
N_vett = ones(1,JJ1)*Inf ;
h_vett = ones(1,JJ1)*Inf ;

% Evaluation of parameters for inversion of LT in each admissible region
find_region = 0 ;
while ~find_region
    for j1 = admissible_regions
        if j1 < J1
            [muj,hj,Nj] = OptimalParam_RB ...
                (t,phi_s_star(j1),phi_s_star(j1+1),p(j1),q(j1),log_epsilon) ;
        else
            [muj,hj,Nj] = OptimalParam_RU(t,phi_s_star(j1),p(j1),log_epsilon) ;
        end
        mu_vett(j1) = muj ; h_vett(j1) = hj ; N_vett(j1) = Nj ;
    end
    if min(N_vett) > 200
        log_epsilon = log_epsilon +log(10) ;
    else
        find_region = 1 ;
    end
end
    

% Selection of the admissible region for integration which involves the
% minimum number of nodes 
[N, iN] = min(N_vett) ; mu = mu_vett(iN) ; h = h_vett(iN) ;

% Evaluation of the inverse Laplace transform
k = -N : N ;
u = h*k ;
z = mu*(1i*u+1).^2 ;
zd = -2*mu*u + 2*mu*1i ;
zexp = exp(z*t) ;
F = z.^(alpha*gama-beta)./(z.^alpha - lambda).^gama.*zd ;
S = zexp.*F ;
Integral = h*sum(S)/2/pi/1i ;

% Evaluation of residues
ss_star = s_star(iN+1:end) ;
Residues = sum(1/alpha*(ss_star).^(1-beta).*exp(t*ss_star)) ;

% Evaluation of the ML function
E = Integral + Residues ;
if isreal(lambda) 
    E = real(E) ;
end

end


% =========================================================================
% Finding optimal parameters in a right-bounded region
% =========================================================================
function [muj,hj,Nj] = OptimalParam_RB ...
    (t, phi_s_star_j, phi_s_star_j1, pj, qj, log_epsilon)

% Definition of some constants
log_eps = -36.043653389117154 ; % log(eps)
fac = 1.01 ;
conservative_error_analysis = 0 ;

% Maximum value of fbar as the ration between tolerance and round-off unit
f_max = exp(log_epsilon - log_eps) ;

% Evaluation of the starting values for sq_phi_star_j and sq_phi_star_j1
sq_phi_star_j = sqrt(phi_s_star_j) ;
threshold = 2*sqrt((log_epsilon - log_eps)/t) ;
sq_phi_star_j1 = min(sqrt(phi_s_star_j1), threshold - sq_phi_star_j) ;

% Zero or negative values of pj and qj
if pj < 1.0e-14 && qj < 1.0e-14
    sq_phibar_star_j = sq_phi_star_j ;
    sq_phibar_star_j1 = sq_phi_star_j1 ;
    adm_region = 1 ;
end

% Zero or negative values of just pj
if pj < 1.0e-14 && qj >= 1.0e-14
    sq_phibar_star_j = sq_phi_star_j ;
    if sq_phi_star_j > 0
        f_min = fac*(sq_phi_star_j/(sq_phi_star_j1-sq_phi_star_j))^qj ;
    else
        f_min = fac ;
    end
    if f_min < f_max
        f_bar = f_min + f_min/f_max*(f_max-f_min) ;
        fq = f_bar^(-1/qj) ;
        sq_phibar_star_j1 = (2*sq_phi_star_j1-fq*sq_phi_star_j)/(2+fq) ;
        adm_region = 1 ;
    else
        adm_region = 0 ;
    end
end

% Zero or negative values of just qj
if pj >= 1.0e-14 && qj < 1.0e-14
    sq_phibar_star_j1 = sq_phi_star_j1 ;
    f_min = fac*(sq_phi_star_j1/(sq_phi_star_j1-sq_phi_star_j))^pj ;
    if f_min < f_max
        f_bar = f_min + f_min/f_max*(f_max-f_min) ;
        fp = f_bar^(-1/pj) ;
        sq_phibar_star_j = (2*sq_phi_star_j+fp*sq_phi_star_j1)/(2-fp) ;
        adm_region = 1 ;
    else
        adm_region = 0 ;
    end
end

% Positive values of both pj and qj
if pj >= 1.0e-14 && qj >= 1.0e-14
    f_min = fac*(sq_phi_star_j+sq_phi_star_j1)/...
        (sq_phi_star_j1-sq_phi_star_j)^max(pj,qj) ;
    if f_min < f_max
        f_min = max(f_min,1.5) ;
        f_bar = f_min + f_min/f_max*(f_max-f_min) ;
        fp = f_bar^(-1/pj) ;
        fq = f_bar^(-1/qj) ;
        if ~conservative_error_analysis
            w = -phi_s_star_j1*t/log_epsilon ;
        else
            w = -2*phi_s_star_j1*t/(log_epsilon-phi_s_star_j1*t) ;
        end
        den = 2+w - (1+w)*fp + fq ;
        sq_phibar_star_j = ((2+w+fq)*sq_phi_star_j + fp*sq_phi_star_j1)/den ;
        sq_phibar_star_j1 = (-(1+w)*fq*sq_phi_star_j ...
            + (2+w-(1+w)*fp)*sq_phi_star_j1)/den ;
        adm_region = 1 ;
    else
        adm_region = 0 ;
    end
end

if adm_region
    log_epsilon = log_epsilon  - log(f_bar) ;
    if ~conservative_error_analysis
        w = -sq_phibar_star_j1^2*t/log_epsilon ;
    else
        w = -2*sq_phibar_star_j1^2*t/(log_epsilon-sq_phibar_star_j1^2*t) ;
    end
    muj = (((1+w)*sq_phibar_star_j + sq_phibar_star_j1)/(2+w))^2 ;
    hj = -2*pi/log_epsilon*(sq_phibar_star_j1-sq_phibar_star_j)...
        /((1+w)*sq_phibar_star_j + sq_phibar_star_j1) ;
    Nj = ceil(sqrt(1-log_epsilon/t/muj)/hj) ;
else
    muj = 0 ; hj = 0 ; Nj = +Inf ;
end

end

% =========================================================================
% Finding optimal parameters in a right-unbounded region
% =========================================================================
function [muj,hj,Nj] = OptimalParam_RU (t, phi_s_star_j, pj, log_epsilon)

% Evaluation of the starting values for sq_phi_star_j
sq_phi_s_star_j = sqrt(phi_s_star_j) ;
if phi_s_star_j > 0
    phibar_star_j = phi_s_star_j*1.01 ;
else
    phibar_star_j = 0.01 ;
end
sq_phibar_star_j = sqrt(phibar_star_j) ;

% Definition of some constants
f_min = 1 ; f_max = 10 ; f_tar = 5 ;

% Iterative process to look for fbar in [f_min,f_max]
stop = 0 ;
while ~stop
    phi_t = phibar_star_j*t ; log_eps_phi_t = log_epsilon/phi_t ;
    Nj = ceil(phi_t/pi*(1 - 3*log_eps_phi_t/2 + sqrt(1-2*log_eps_phi_t))) ;
    A = pi*Nj/phi_t ;
    sq_muj = sq_phibar_star_j*abs(4-A)/abs(7-sqrt(1+12*A)) ;
    fbar = ((sq_phibar_star_j-sq_phi_s_star_j)/sq_muj)^(-pj) ;
    stop = (pj < 1.0e-14) || (f_min < fbar && fbar < f_max) ;
    if ~stop
        sq_phibar_star_j = f_tar^(-1/pj)*sq_muj + sq_phi_s_star_j ;
        phibar_star_j = sq_phibar_star_j^2 ;
    end
end
muj = sq_muj^2 ;
hj = (-3*A - 2 + 2*sqrt(1+12*A))/(4-A)/Nj ;

% Adjusting integration parameters to keep round-off errors under control
log_eps = log(eps) ; threshold = (log_epsilon - log_eps)/t ;
if muj > threshold
    if abs(pj) < 1.0e-14 , Q = 0 ; else Q = f_tar^(-1/pj)*sqrt(muj) ; end
    phibar_star_j = (Q + sqrt(phi_s_star_j))^2 ;
    if phibar_star_j < threshold
        w = sqrt(log_eps/(log_eps-log_epsilon)) ;
        u = sqrt(-phibar_star_j*t/log_eps) ;
        muj = threshold ;
        Nj = ceil(w*log_epsilon/2/pi/(u*w-1)) ;
        hj = sqrt(log_eps/(log_eps - log_epsilon))/Nj ;
    else
        Nj = +Inf ; hj = 0 ;
    end
end

end