In boron neutron capture therapy (BNCT)-a form of binary radiotherapy-the primary challenge in treatment planning systems for dose calculations arises from the time-consuming nature of the Monte Carlo (MC) method. Recent progress, including the use of neural networks (NN), has been made to accelerate BNCT dose calculations. However, this approach may result in significant dose errors in both the tumor and the skin, with the latter being a critical organ in BNCT. Furthermore, owing to the lack of physical processes in purely NN-based approaches, their reliability for clinical dose calculations in BNCT is questionable. In this study, a physically constrained MC-NN (PCMC-NN) coupling algorithm is proposed to achieve fast and accurate computation of the BNCT three-dimensional (3D) therapeutic dose distribution. This approach synergizes the high precision of the MC method with the speed of the NN and utilizes physical conservation laws to constrain the coupling process. It addresses the time-consuming issue of the traditional MC method while reducing dose errors. Clinical data were collected from 113 glioblastoma patients. For each patient, the 3D dose distributions for both the coarse and detailed dose grids were calculated using the MC code PHITS. Among these patients, the data from 14 patients were allocated to the test set, 9 to the validation set, and the remaining to the training set. A neural network, 3D-Unet, was built based on the coarse grid dose and patient CT information to enable fast and accurate computation of the 3D detailed grid dose distribution of BNCT. Statistical evaluations, including relative deviation, dose deviation, mean absolute error (MAE), and mean absolute percentage error (MAPE) were conducted. Our findings suggested that the PCMC-NN algorithm substantially outperformed the traditional NN and interpolation methods. Furthermore, the proposed algorithm significantly reduced errors, particularly in the skin and GTV, and improved computational accuracy (hereinafter referred to simply as 'accuracy') with a MAPE range of 1.6%-4.0% and a maximum MAE of 0.3Gy (IsoE) for different organs. The dose-volume histograms generated by the PCMC-NN aligned well with those obtained from the MC method, further validating its accuracy. The PCMC-NN algorithm enhanced the speed and accuracy of BNCT dose calculations by combining the MC method with the NN algorithm. This indicates the significant potential of the proposed algorithm for clinical applications in optimizing treatment planning.