This paper is intended to present an implementation of the hypersingular boundary integral equations in terms of first-order stress functions for stress computations in plane orthotropic elasticity. In general, the traditional computational technique of the boundary element method used for computing the stress distribution on the boundary and close to it is not as accurate as it should be. In contrast, the accuracy of stress computations on the boundary is greatly increased by applying the hypersingular integral equations. Contrary to the method in which the solution is based on an approximation of displacement field, here the first-order stress functions and the rigid body rotation are the fundamental variables. An advantage of this approach is that the stress components can be obtained directly from the stress functions, there is, therefore, no need for Hooke's law, which should be used when they are computed from displacements. In addition, the computational work can be reduced when the stress distribution is computed at an arbitrary point on the boundary. The numerical examples presented prove the efficiency of this technique.