We have developed a method and algorithms for solving the generalized Lyapunov equation for a wide class of continuous time-varying bilinear systems based on the Gramian method and an iterative solution construction method proposed earlier for such equations. The approach consists in diagonalizing the original system, obtaining a separable spectral decomposition of the Gramian of the time-invariant linear part in terms of the combination spectrum of the dynamics matrix of the linear part, applying the spectral decomposition of the kernel matrix of the solution obtained at the previous step at each iteration step, and then aggregating the matrix entries. A spectral decomposition of the Gramians of controllability and observability of a time-varying bilinear system is obtained as the sum of sub-Gramian matrices corresponding to pair combinations of the eigenvalues of the dynamics matrix of the linear part. A new method and algorithm for entry-by-entry calculation of matrices for solving the generalized Lyapunov equation for bilinear systems has been developed. The fundamental novelty of the approach lies in the transfer of calculations from the solution matrix to the calculation of the sequence of its entries at each iteration step.