A higher order semi-analytical method for the simulation of heat transfer in fields of geothermal boreholes is presented. The method uses the integration of the point heat source solution to evaluate borehole wall and ground temperatures. Instead of piecewise uniform heat extraction rate and borehole wall temperature along segments of discretized boreholes, the presented method considers polynomial variations of these quantities using the superposition of polynomial basis functions. Increasing the order of the polynomial basis function (by increasing the number of nodes per borehole segment) is shown to have a greater impact on the model accuracy than increasing the number of segments. The error on the g-function of a field of 8 inclined boreholes is less than 0.1 % using 2 segments of 11 nodes. The presented method enables the simulation of boreholes with curved trajectories. The results show an error of 10.3 % between the thermal response evaluated using the real measured trajectories of boreholes in a real installation and the thermal response evaluated using approximated straight trajectories.