Coalescent methods are proven and powerful tools for population genetics, phylogenetics, epidemiology, and other fields. The multispecies coalescent (MSC) model has been widely employed by phylogenetic algorithms to construct the species tree while accounting for incomplete lineage sorting (ILS). However, the no-recombination assumption of the MSC model has been questioned. To analyze large genomic regions, we need to simultaneously account for both ILS and recombination. A promising avenue for the analysis of large genomic alignments, which are now commonplace, are coalescent hidden Markov model (coalHMM) methods, but these methods have lacked general usability and flexibility. I introduce in this thesis a novel method, VICAR (Variational Inference under the CoAlescent with Recombination), for automatically learning a coalHMM and inferring the posterior distributions of evolutionary parameters using black-box variational inference, with the transition rates between local genealogies derived empirically by simulation. This derivation enables VICAR to work directly with three or four taxa and through a divide-and-conquer approach with more taxa. Using a simulated data set resembling a human-chimp-gorilla scenario, I show that VICAR has comparable or better accuracy to previous coalHMM methods. Both species divergence times and population sizes were accurately inferred. The method also infers local genealogies and I report on their accuracy. Furthermore, I illustrate how to scale the method to larger data sets through a divide-and-conquer approach. This accuracy means that my approach is useful now, and by deriving transition rates by simulation it is flexible enough to enable future implementations of all kinds of population models. I have implemented VICAR in the publicly available software package PhyloNet.