An algorithm plotting contour lines for discrete values z ij , given at the nodes of a rectangular mesh is described. A bicubic Hermite polynomial f ( x , y ) is determined for every rectangle of the mesh, interpolating the z ij and the derivatives z x , z y , and z xy . The derivatives are optionally computed by the algorithm. The contours found are normally smooth curves. They consist of polygons approximating intersections with the bicubics. It is possible to fill the areas between them with certain colors or patterns. This is done with a piecewise technique rectangle by rectangle. The method for finding the points of the polygons is shortly reviewed, and some numerical problems are pointed out. The algorithm has a flexible, easy-to-use interface and is easily installed with all plotting systems, provided that a fill-area command is available. A GKS interface may be used.