This paper proposes a global iterative sliding mode control approach for high-precision contouring tasks of a flexurelinked biaxial gantry system. For such high-precision contouring tasks, it is the typical situation that the involved multi-axis cooperation is one of the most challenging problems. As also would be inevitably encountered, various factors render the multi-axis cooperation rather difficult; such as the strong coupling (which naturally brings nonlinearity) between different axes due to its mechanical structure, the backlash and deadzone caused by the friction, and the difficulties in system identification, etc. To overcome the above-mentioned issues, this work investigates an intelligent model-free contouring control method for such a multiaxis motion stage. Essentially in the methodology developed here, it is firstly ensured that all the coupling, friction, nonlinearity, and disturbance (regarded as uncertain dynamics in each axis) are suitably posed as 'uncertainties'. Then, a varying-gain sliding mode control method is proposed to adaptively compensate for the matched unknown dynamics in the time domain, while an iterative learning law is applied to suppress the undesirable effects (arising from the repetitive matched and unmatched uncertainties in the iteration domain). With this approach, the chattering that typically results from the overestimated control gains in the sliding mode control is thus suppressed during the iterations. To analyze the contouring performance and show the improved outcomes, rigorous proof is furnished on both the stability in the time domain and the convergence in the iteration domain; and the real-time experiments also illustrate that the requirements of precision motion control towards highspeed and complex-curvature references can be satisfied using the proposed method, without prior knowledge of the boundary to the unknown dynamics.