Subdivision surface is an efficient geometric modeling tool that is generated by repeated approximation or interpolation from an initial control mesh. For each time of subdivision the old mesh will be split and the vertexes of the mesh will be refined for a new mesh. Most existing algorithms aim to construct a smooth subdivision surface by smoothing the vertexes of the split mesh. When some features are desired, special rules or parameters should often be applied. In this paper a new non-stationary subdivision method is presented for surface modeling by recursive interpolation from an initial triangular control mesh. For each time of subdivision the mesh itself is a piecewise linear approximation to the final surface, and more details will be added to the surface when the triangles have been split and replaced by more refined sub-triangles. We compute every interpolated edge vertex as a solution to the problem of least square fitting of the edge and its neighboring triangles. Sharp features as well as flat regions can be kept well during subdivision. To obtain a visually smooth subdivision surface in the limit, the normal at the vertexes of old meshes will be used as constraint for every time of subdivision. The main advantage of this new method is that the interpolated vertexes depend on the local geometry of the mesh, but not the valances of the mesh. The examples we have tested show that the interpolating surfaces by this new method can inherit the shape of initial course mesh more fairly and naturally.