    for jj=1:4
        % Caclulate time derivative of states
        xx = x(:, jj);
        f = [
            xx(2)
            -omega_n^2*xx(1) - 2*zetaa*omega_n*xx(2) + omega_n^2*u
            ];

        % Gill's variables
        k(:, jj) = a(jj) * (f - b(jj)*q(:, jj));
        x(:, jj+1) = x(:, jj) + k(:, jj)*dt;
        q(:, jj+1) = q(:, jj) + 3*k(:, jj) - c(jj)*f;
    end
    % Update Gill's variables for use in the next time step
    q(:, 1) = q(:, 5);
    x(:, 1) = x(:, 5);