The homotopy method is implemented as the first applicable method for directly solving three dimensional frictional contact problems with the orthotropic Coulomb friction. A system of nonlinear equations are derived from the nonlinear complementary relations which constitute an exact description of the three-dimensional frictional contact phenomenon. No linearlization of the elliptic friction law, previously introduced to derive approximated linear complemetarity problem, is introduced. The equations are presented in an incremental form using updated Lagrangian description to accommodate the loading-path dependency of friction phenomena and large displacement cases. The probability-one homotopy method known as a globally convergent zero-finding algorithm is then adopted for a solution scheme of this formulation and applied to each incremental step. An ordinary differential equation based method is used as a zero curve following method. Several three dimensional frictional contact problems are solved to test the method and the results compared with those of commercial package and other approximations. The first one is a two body contact between an elastic punch and a rigid half surface. The results are compared with those of ABAQUS, a general purpose nonlinear analysis package, where a penalty method has been used. Good agreements are seen for the normal pressure distributions and tangential displacements. The second example is a contact between a rigid punch and an elastic half surface handled by other author using octagonal approximation of the friction law. Good agreements in the normal pressure distributions during the loading and unloading are observed. But, the directions and the amount of tangential displacements are observed somewhat different, which is thought to be caused by the octagonal approximation of the friction law in the literature. The subsequent development of outward slip-stick-inward slip motion during unloading process mentioned by several authors for two dimensional and axisymmetric cases are well observed. The last example is a three body contact between two elastic bodies and a rigid surface and rigid body motions are considered for each elastic body. Unlike the above two examples, changes in contact area during loading and unloading process occur in this example and have been well simulated. The detail phenomeon of stick-slip motions and pressure distributions on each contact area are shown relatively well. The results at various friction coefficients including orthotropic cases are also presented. In summary, the method is shown to be a viable strong tool for detailed and precise analysis of complicated phenomena which involve three dimensional frictional contact. The numerical efficiency is much dependent on the zero curve following algorithm which needs futher attention.