Does using Multifrequency eliminate totally the ionosphere delay on GPS measurements?
As @CamilleGoudeseune warns, no nothing is perfect except mathematics, though I'd add LeeLoo to the list.
With only two frequencies the correction is pretty good but not perfect.
Dispersion is everywhere. Whenever there is the possibility of absorption at one wavelength, we see a frequency dependent phase shift at many wavelengths.
In this case I believe that it is absorption at low frequencies (below the MUF or Maximum Usable Frequency (for ionospheric reflection from below) at say 30 to 60 MHz) that results in frequency-dependent delays (dispersion) in the low GHz region, but precipitable water vapor is also a player and this is the reason that radio astronomers also use wavefront correction for large interferometric arrays.
Wikipedia's Global Positioning System; Satellite frequencies says:
All satellites broadcast at the same two frequencies, 1.57542 GHz (L1 signal) and 1.2276 GHz (L2 signal).
and
The L4 band at 1.379913 GHz is being studied for additional ionospheric correction.
Radio astronomers use dispersion from the interstellar medium to estimate the distance at which fast radio burst (FRB) events happen. In this case the signal is broadband and you can plot exactly the arrival time versus frequency and fit a pretty good dispersion model to it.
Currently using only the two L1 and L2 frequencies, a GPS receiver algorithm can give only a "first order approximation" to ionospheric and other dispersion effects. The L4 frequency (mentioned above) between the two might potentially be used to improve the correction.
For further reading see Impact and Implementation of Higher-Order Ionospheric Effects on Precise GNSS Applications
Waterfall plot illustration of how for a given signal dispersion results in increasingly delayed arrival times for lower frequencies, as they get closer to the plasma frequency. From the linked answer.
From Lorimer et al. (cited above):

Figure 2: Frequency evolution and integrated pulse shape of the radio burst. The survey data, collected on 2001 August 24, are shown here as a two-dimensional ‘waterfall plot’ of intensity as a function of radio frequency versus time. The dispersion is clearly seen as a quadratic sweep across the frequency band, with broadening towards lower frequencies. From a measurement of the pulse delay across the receiver band using standard pulsar timing techniques, we determine the DM to be 375±1 cm−3 pc. The two white lines separated by 15 ms that bound the pulse show the expected behavior for the cold-plasma dispersion law assuming a DM of 375 cm−3 pc. The horizontal line at ∼ 1.34 GHz is an artifact in the data caused by a malfunctioning frequency channel. This plot is for one of the offset beams in which the digitizers were not saturated. By splitting the data into four frequency sub-bands we have measured both the half-power pulse width and flux density spectrum over the observing bandwidth. Accounting for pulse broadening due to known instrumental effects, we determine a frequency scaling relationship for the observed width W = 4.6 ms (f/1.4 GHz)−4.8±0.4 , where f is the observing frequency.
A power-law fit to the mean flux densities obtained in each sub-band yields a spectral index of −4 ± 1. Inset: the total-power signal after a dispersive delay correction assuming a DM of 375 cm−3 pc and a reference frequency of 1.5165 GHz. The time axis on the inner figure also spans the range 0–500 ms.