<p>Ionospheric function fitting model can play a pivotal role in addressing the issue of understanding ionospheric geographical variations. However, the optimal functional fitting form of an ionospheric model has not yet been convincingly determined so far. Since the parameters to drive NeQuick model are available in broadcast ephemeris and NeQuick model has better accuracy than other broadcast ionospheric models, NeQuick model was employed to determine the structure of ionospheric function fitting model by Akaike information criterion in advance, in which the order of polynomials along latitude and longitude was determined from the candidate sets for each period. Global navigation satellite system data from continuously operating reference stations in Hong Kong and the United States, collected on day of the year (DOY) 102, 2021 and DOY 225, 2024, were utilized to establish the flexible model. To validate the model, model accuracy and differential slant total electron content (dSTEC) observations at training and checking stations were analyzed. Compared with the traditional polynomial function, the accuracy of the flexible model in terms of training errors, fitting errors, and dSTEC observations has been improved by up to 40%, 26%, and 46%, respectively, under disturbed ionospheric conditions. In addition, the flexible model constructed for the calmer ionospheric conditions gain a slight accuracy improvement in all three aspects. The findings gathered together indicate that the flexible model is capable of capturing the ionospheric finer variations.</p>