A polarizable model potential (PMP) function for methanol and three ions (Cl-, Na+, and Mg2+) is developed on the basis of high-level ab initio molecular orbital calculations (MP2/6-31++G(2d,2p)). The PMP function consists of Coulomb, polarization, and Lennard-Jones terms. The permanent partial charge parameters of the Coulomb term, and the multicenter polarizability parameters of the polarization term, are determined by using electrostatic potential optimization, and polarizable one-electron potential (POP) optimization, respectively. The Lennard-Jones parameters an adjusted to reproduce the ab initio potential energy surfaces of methanol dimer and ion-methanol complexes. The PMP function using the final parameters reproduced well the ab initio energy surfaces of methanol dimer, ion-methanol, and methanol-ion-methanol systems. The root-mean-square (mms) deviations of ion-methanol systems having various conformations are only 10-12% of the ab initio interaction energies. The electron density changes by the complex formations were reproduced well. The rms deviations are 1-2% of the surface electrostatic potentials. The interaction energies of each classical term were compared with the corresponding terms derived from an ab initio energy decomposition analysis (EDA). The Coulomb and polarization energy reproduced well the electrostatic and induction energy of EDA, respectively. The Lennard-Jones energy reproduced well the sum of EDA exchange repulsion energy, dispersion energy, and deformation energy. In the PMP function, most many-body effects can be adequately evaluated by the introduction of multicenter polarizabilities determined from the POP optimization.