In this paper, an efficient computational method is presented for state space analysis of bilinear systems via Legendre wavelets. The differential matrix equation is converted to a generalized Sylvester matrix equation by using Legendre wavelets as a basis. First, an explicit expression for the inverse of the integral operational matrix of the Legendre wavelets is presented. Then using it, we propose a preorder traversal algorithm to solve the generalized Sylvester matrix equation, which greatly reduces the computation time. Finally the efficiency of the proposed method is discussed using numerical examples.