    function [t, y, Z] = fcoll(f, b, gam, alpha, eta, r, N)
%
%  [t, y] = fcoll(f, b, gam, alpha, eta, r, N)%
%
%  Code from:
%
%  Cardone, A., Conte, D., Paternoster, B. 
%  A MATLAB implementation of spline collocation methods for 
%  fractional differential equations. 
%  Lect. Notes Comput. Sci. 12949 LNCS, 387-401 (2021).
%

    % Time grid with graded mesh
    t = b * ((0:N)' / N).^r;
    h = diff(t);
    
    % Initialize matrices and vectors
    Alagr = matrix_Lagrange(eta);
    A = matrix_A(alpha, eta, Alagr);
    bvect = vector_b(alpha, Alagr);
    
    % Solver options
    options = optimset('TolFun', 1e-14, 'TolX', 1e-14, 'Display', 'off');
    
    m = length(eta);
    Z = zeros(m, N);
    y = zeros(N+1, 1);
    y(1) = gam(1);
    
    % Initial guess
    iniz = feval(f, t(1) + eta*h(1), gam(1)*ones(m, 1));
    
    % Time-stepping loop
    for j = 1:N
        tj = t(j) + eta*h(j);
        
        % Compute history terms
        Bj = lag(Z, Alagr, t, h, eta, alpha, j);
        Qj = Q(tj, alpha, gam);
        
        % Solve nonlinear system
        Z(:, j) = fsolve(@(x) system_F(x, f, tj, A, Bj, Qj, h, j, alpha), iniz, options);
        
        % Update solution
        Wj = lag_y(Z, Alagr, t, h, alpha, j);
        Qjp1 = Q(t(j+1), alpha, gam);
        y(j+1) = h(j)^alpha * (bvect') * Z(:, j) + Wj + Qjp1;
        
        % Update initial guess for next step
        iniz = Z(:, j);
    end
end