Memristors have emerged as a promising candidate for next-generation nonvolatile memory and neuromorphic computing architectures due to their inherent capability of multilevel resistive switching and non-linear current-voltage (I–V) behavior. A crucial challenge in advancing their practical utility lies in accurately modeling their electrical characteristics, especially the pinched hysteresis loops (PHLs) observed during bipolar voltage sweeps. While prior models have successfully captured single PHLs, the simulation of double PHLs, essential for representing multilevel resistance states, remains a relatively unexplored domain. This chapter details the development, implementation, and optimization of a Verilog-A-based simulation framework for modeling double-PHL behavior in Au/TiO2/Ti-based resistive random-access memory (RRAM) devices. Through a step-by-step methodology, incorporating the Simmons tunneling model, we introduce parameter modifications, dual-threshold voltage insertion, and segmented simulation regimes to reflect experimental I–V behavior with improved accuracy. The optimized simulation results demonstrate high fidelity to experimental data, with error margins maintained within 10%. This chapter provides a crucial simulation foundation for implementing multilevel RRAM devices in real-time circuits.