Approximating the partition function of the ferromagnetic Ising model with general external fields is known to be #BIS-hard in the worst case, even for bounded-degree graphs, and it is widely believed that no polynomial-time approximation scheme exists. This motivates an average-case question: are there classes of instances for which polynomial-time approximation schemes exist? We investigate this question for the random field Ising model on graphs with maximum degree $Δ$. We establish the existence of fully polynomial-time approximation schemes and samplers with high probability over the random fields if the external fields are IID Gaussians with variance larger than a constant depending only on the inverse temperature and $Δ$. The main challenge comes from the positive density of vertices at which the external field is small. These regions, which may have connected components of size $Θ(\log n)$, are a barrier to algorithms based on establishing a zero-free region, and cause worst-case analyses of Glauber dynamics to fail. The analysis of our algorithm is based on percolation on a self-avoiding walk tree.