This paper outlines a numerical algorithm capable of simulating the full three-dimensional dynamics of magnetic fluid droplets in external magnetic fields by solving boundary integral equations. The algorithm works with arbitrary droplet and carrier fluid viscosity ratios. It is validated with known theoretical relationships. It also enables evaluating various approximations often used to describe ellipsoidal droplets by comparing the droplet dynamics calculated from such approximations to the results obtained from first principles using our numerical algorithm. The algorithm may be used to investigate droplet configurations in arbitrary magnetic fields and to indirectly calculate the physical properties of magnetic fluid droplets and predicting the magnetic field thresholds above which the droplet shape can develop instabilities in the form of various spikes.