An improved response matrix method has been proposed to effectively take into account the anisotropy of neutron angular distributions. The method utilizes a relation between the P0 and P1 components of a neutron angular distribution instead of calculating them independently. Hence the number of unknowns as well as computing time can be kept about the same as in the conventional response matrix method which adopts an isotropic approximation of a neutron angular distribution. The proposed method has been evaluated by applying it to one-dimensional slab and two-dimensional hexagonal systems. The results are quite promising: In comparison with the reference SN calculation, the difference of the neutron multiplication factor and power distribution is within 0.1% Δk/k and 2%, respectively, and furthermore, the computing time is reduced to below one-third.