In 1998, Shibata and Sasaki [Phys. Rev. D 58, 104011 (1998)] presented an approximate analytical formula for the radius of the innermost stable circular orbit (ISCO) of a neutral test particle around a massive, rotating and deformed source. In the present paper, we generalize their expression by including the magnetic dipole moment. We show that our approximate analytical formulas are accurate enough by comparing them with the six-parametric exact solution calculated by Pach\'on et. al. [Phys. Rev. D 73, 104038 (2006)] along with the numerical data presented by Berti and Stergioulas [MNRAS 350, 1416 (2004)] for realistic neutron stars. As a main result, we find that in general, the radius at ISCO exhibits a decreasing behavior with increasing magnetic field. However, for magnetic fields below 100GT the variation of the radius at ISCO is negligible and hence the non-magnetized approximate expression can be used. In addition, we derive approximate analytical formulas for angular velocity, energy and angular momentum of the test particle at ISCO.