Stochastic Differential Equations (SDEs) are used as statistical models in many disciplines. However, intractable likelihood functions for SDEs make inference challenging, and we need to resort to simulation-based techniques to estimate and maximize the likelihood function. While importance sampling methods have allowed for the accurate evaluation of likelihoods at fixed parameter values, there is still a question of how to find the maximum likelihood estimate. We propose an efficient Gaussian-process-based method for exploring the parameter space using estimates of the likelihood from an importance sampler. Our technique accounts for the inherent Monte Carlo variability of the estimated likelihood, and does not require knowledge of gradients. The procedure adds potential parameter values by maximizing the so-called expected improvement, leveraging the fact that the likelihood function is assumed to be smooth. Our simulations demonstrate that our method has significant computational and efficiency gains over existing grid- and gradient-based techniques. Our method is applied to the estimation of ocean circulation from Lagrangian drift data in the South Atlantic ocean.