A two-dimensional inverse method is presented for estimating the boundary heat flux imposed on phase change materials in the melting/solidification process. The inverse problem is handled by the conjugate gradient method with the adjoint problem for function estimation. The method is verified using unsynthetic data, and the effect of the number of sensors, sensor placement, and measurement errors on the accuracy is investigated. In all cases, it is found that the method remains stable and accurate with a root mean square relative error less than 10%. Also, the number of sensors depends on the desired accuracy, computing time, and frequency of spatial changes in heat flux. Using field measurement is the best option for high accuracy and time-efficient estimation of spatial changes in the heat flux function. The accuracy of the estimated heat flux strongly depends on the sensor distance from the boundary. Therefore, it is desirable to place the sensor as close to the boundary as possible. With the addition of measurement errors, the presented method remains stable even with high measurement errors. This shows that this method can be used for real conditions.