function [y1,T,y0] = prob3(t,y,jaco)
%
% [alpha,T,y0] = prob3
%
% or
%            f = prob3(t,y)
%
% or
%           f' = prob3(t,y,1)
%
% or
%            y = prob3(t)
%
% Problem 3 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;
A    = [-92 -87; -58 -63]/5; 
V    = [1 1.5; -1 1];  % A = V*diag([-1 -30])/V
b    = [67 83]/10;
if nargin==0    % problem data
    y1   = alfa;
    T    = 100;
    y0   = [5 10];
elseif nargin==1 % solution
    V1      = inv(V)';
    ybar    = A\b';
    ybar    = ybar';
    y0b     = ([5 10]-ybar)*V1;
    talf    = t(:).^alfa;
    ml1     = ml(-talf,alfa);
    ml2     = ml(-30*talf,alfa);
    e       = ones(length(t),1);
    y1      = kron(e,y0b);
    y1(:,1) = y1(:,1).*ml1;
    y1(:,2) = y1(:,2).*ml2;
    y1      = y1*V' + ybar; 
elseif nargin==2 % vector field
    if length(t)==1
       y1 = A*y(:) -b';
    else
       y1 = y*A' -kron(ones(length(t),1),b);
    end   
else             % Jacobian
    y1    = A;
end
return
