function [ut] = ComputeControl(ct, s_int, kd)
    % Integrate controls using the Trapezoidal rule and approxmiate derivatives.
    % kd is a function handle to a function that calculates kd - kdsquared
    % for delta^2 CCM, kdfourth for delta^4 CCM.
    
    ut = zeros(1/s_int + 1,1);
    for s = s_int:s_int:1
        e = round(s/s_int + 1); % end index for current s
        % calculate end derivatives using skewed stencils
        derivative0 = (ct(2) - ct(1))/s_int;
        derivativen = (ct(e) - ct(e - 1))/s_int;
        ut(e) = 0;
        ut(e) = ut(e) + kd(derivative0);
        ut(e) = ut(e) + kd(derivativen);
        for i = 2:e-1
            % calculate derivative using centered stencil
            derivative = (ct(i + 1) - ct(i - 1))/2/s_int;
            ut(e) = ut(e) + 2*kd(derivative);
        end
        ut(e) = s_int/2 * ut(e);
    end
end