Drift-diffusion models that account for the motion of both electronic and ionic charges are important tools for explaining the hysteretic behaviour. Furnishing numerical solutions to such models for realistic operating conditions is challenging owing to the extreme values of some of the parameters. We present a finite difference scheme with time step that provides second order accuracy in the mesh spacing. The method is able to use realistic parameters values whilst providing high accuracy. Also, three diagonal matrix approximation (TDMA) method is exploit due to matrix solve simplification. This method is robust and fast way to solve. Ion vacancy density, electron concentration, and hole concentration profiles are calculated in transient time-scale. In following, built-in potential is varied and profiles are illustrated. This approach paves the way to have a better insight of device physics and its related phenomena such as ionic motion, hysteresis.