function [y1,T,y0] = prob4(t,y,jaco)
%
% [alpha,T,y0] = prob4
%
% or
%            f = prob4(t,y)
%
% or
%           f' = prob4(t,y,1)
%
% or
%            y = prob4(t)
%
% 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