function [t,y] = tsfcoll(f, b, gam, alpha, eta, r, N)
%
%  [t,y] = tsfcoll(f, b, gam, alpha, eta, r, N)
%
%  Code from:
%
%  A.Cardone, D.Conte, B.Paternoster.
%  A Matlab code for fractional differential equations based on two-step
%  spline collocation methods.
%  In "Fractional Differential Equations", Springer, 2023, pp. 121-146.
%
%
    t = b * ((0:N)' / N).^r;
    h = diff(t);
    m = length(eta);
    etatilde = linspace(0, 1, 2*m)';
    
    % Starting solution in [t(1), t(2)]
    [Alagr_tilde] = matrix_Lagrange(etatilde);
    [~, ytilde, ztilde] = fcoll(f, t(2), gam, alpha, etatilde, r, 1);
    
    H = [fliplr(vander(eta)), diag(eta.^m) * fliplr(vander(eta))];
    Z = zeros(m, N); 
    Z(:, 1) = H * Alagr_tilde' * ztilde;
    y = zeros(N+1, 1); 
    y(1) = gam(1); 
    y(2) = ytilde(2);
    
    Alagr = zeros(2*m, 2*m, N);
    options = optimset('TolFun', 1e-14, 'TolX', 1e-14, 'Display', 'off');
    
    for j = 2:N
        tj = t(j) + eta * h(j);
        Alagr(:, :, j) = tsmatrix_Lagrange(eta, h, j);
        A = tsmatrix_A(alpha, eta, Alagr(:, :, j));
        E = tsmatrix_E(alpha, eta, Alagr, t, h, j);
        Bj = tslag(Z, ztilde, A, E, Alagr, Alagr_tilde, t, h, eta, alpha, j);
        Qj = Q(tj, alpha, gam);
        iniz = Z(:, j-1);
        Z(:, j) = fsolve(@tssystem_F, iniz, options, f, tj, A, Bj, Qj, h, j, alpha);
        
        b = tsmatrix_A(alpha, 1, Alagr(:, :, j));
        g = tsmatrix_E(alpha, 1, Alagr, t, h, j);
        Wj = tslag(Z, ztilde, b, g, Alagr, Alagr_tilde, t, h, 1, alpha, j);
        Qjp1 = Q(t(j+1), alpha, gam);
        y(j+1) = h(j)^alpha * b(:, m+1:2*m) * Z(:, j) + Wj + Qjp1;
    end
end