A multiple-image method is proposed to approximate the reaction-field potential of a source charge inside a finite length cylinder due to the electric polarization of the surrounding membrane and bulk water. When applied to a hybrid ion-channel model, this method allows a fast and accurate treatment of the electrostatic interactions of protein with membrane and solvent. To treat the channel/membrane interface boundary conditions of the electric potential, an optimization approach is used to derive image charges by fitting the reaction-field potential expressed in terms of cylindric harmonics. Meanwhile, additional image charges are introduced to satisfy the boundary conditions at the planar membrane interfaces. In the end, we convert the electrostatic interaction problem in a complex inhomogeneous system of ion channel/membrane/water into one in a homogeneous free space embedded with discrete charges (the source charge and image charges). The accuracy of this method is then validated numerically in calculating the solvation self-energy of a point charge.