A new analytical method is presented to compute the seismic response of three-dimensional alluvial valleys. In generalized coordinates along the boundary shape of a three-dimensional alluvial valley, wave field equations and boundary conditions are formulated as a two-point boundary value problem for linear first-order ordinary differential equations. The Riccati method is applied to solve the problem. The numerical results show the validity of the presented method. Complicated amplification patterns in a frequency-space domain induced in alluvial valleys are simply explained by interference phenomena among the direct wave and the surface waves generated on valley edges.