We present a simple method for solving a free boundary value problem of locating either minimum free energy transition path ensemble or minimum potential energy transition path between two configurations on a corresponding energy surface. Our method, called gradient-augmented Harmonic Fourier Beads, employs the global Fourier representation of the path that is a curve interpolation of a discrete set of points on the surface - beads. To optimize the path curve, the method computes energy gradients for each bead from either convex optimization or molecular dynamics or Monte Carlo simulations subject to harmonic restraints. The path optimization is driven by primitive Steepest Descent procedure. Line integration of the Fourier transformed forces along the path curve provides complete and accurate structural and energetic information regarding all the intermediates and saddle points present. The utility of the HFB method is demonstrated by computing potentials of mean force for various transformations of diverse molecular systems.