In order to implement the cross-media gamut mapping in colour image reproduction, the gamut of each medium involved should be first determined. In this paper, a new model was proposed to fit the surface of gamut boundary by two-variable high order polynomials. A Neso cathode ray tube monitor was used to illustrate how to apply the model to determining the medium colour gamut boundary and to evaluate the accuracy of this method. The colour gamut boundary could be analytically expressed by six two-variable high order polynomial equations. The coefficients of each equation could be calculated by least square approximation using the data measured of the colours on the boundary. Several key issues about this method, including the determination accuracy, the optimal quantities of the colour samples and the appropriate orders, were discussed. The experimental data analysis showed that the colour gamut boundary of the CRT could be described accurately using six sets of two-variable high order polynomials, when the order was set to four and 5×5×6 colours on the gamut boundary were used for model training.