We propose a new algorithm to plot principal stress trajectories in plane systems. Much time and labor with experience are required in order to draw trajectories from the result of the photoelastic investigations. The possibility of branching or joining of trajectories cannot be disregarded, when the method of tracking principal stress direction from element to element is used. Our algorithm, on the contrary, plots the principal stress trajectories as the contours of the direction function which is evaluated from the result of FEM using the constant strain triangular elements. Direction function is a two dimensional scalar field whose value is constant along principal stress trajectories. Two function values for the first and the second principal stress of each node are computed by solving the simultaneous liner equations with rectangular sparse matrix. The unit direction vectors of principal stress given by FEM determine the matrix. When the scalar field forms a spiral surface, the unique value of the direction potential of each node is obtained by putting the discontinuity along a polygonal line in the finite element mesh. The line density of the trajectories with broken line pattern in the output figure is read as the magnitude of the principal stress declining toward the free boundary.