Note: The code is run from the main program called RangeExpansion.m Please change the input parameters in RangeExpansion.m Vs and alleffect were varied in "The effect of the recombination rate between adaptive loci on the capacity of a population to expand its range", where: Vs = 2, alleffect = 1/sqrt(20) or, Vs = 0.5, alleffect = 0.2. The user may also want to change: rm (maximal intrinsic growth rate), mu (mutation rate), sig (dispersal distance) or M (the number of demes in the habitat). The parameter b determines the slope of the environmental gradient in the center, please note that once this is set, the slope at the edges and the number of loci are calculated automatically, so that enough loci are required for adaptation to the entire habitat and the expected growth rate at the edges is 0 or smaller. The recombination rate is set using recrate: in the present study it was set to either 0.5, 0.05, 0.01, 0.005, 0.001, 0.0001 or 0.00001. The probability to accept selfing is set using self: in the present study this parameter was either 0 or 1, but it may be set to any value between 0 and 1, meaning that if two gametes are picked at random, the probability to accept them to form a zygote is self. The probability to accept non-selfing gametes is always 1. The function DiploidHabitat(b) automatically creates the environmental gradient, "theta", and its slope "dtheta", using the input above. It also determines the number of loci (parameter called "loci") and the location of the centers of clines ("ccenters"). The parameter g sets the number of generations that are realised. The parameter snap sets the number of generations between measurement points, when output is saved. K1 is the maximum carrying capacity in each deme. In the present study, it was set to 150, but may be changed as desired. InitVar sets the initial standing genetic variation, it can be set either to 0 (meaning that the "low standing genetic variation" simulations are performed) or to 1 (meaning that the "high standing genetic variation" simulations are performed). It cannot be set to other values than 0 or 1. The function InitCond(ccenters,K1,init_range,InitVar) creates the initial conditions required for a realisation, using the parameters created above. The function DiploidMain(mu,recrate,K,N,theta,input_pop,TL) runs the main algorithm, using sub-functions called DiploidSelection, DiploidMutation, DiploidFertilisation and DiploidMigration. DiploidMain generates the following output: output_pop: a cell array explicitly containing each individual with its genotype in each deme at the end of the realisation (i.e. in generation g); please note that this is a logical array, so that each nonzero allele effect is set to 1. N: g/snap by M+1 matrix that contains the local population size in each deme. Npm: the same as N, but with population size saved before migration occurred in the present generation. Ntot: the total population size in the entire habitat at the end of the simulation. Vg: g/snap by M+1 matrix that contains the local LE part of the genetic variance in each deme. Pself: g/snap by M+1 matrix that contains the local realised selfing frequency in each deme. cline: loci by M+1 by 5 matrix containing the allele frequencies at each locus in each deme at 5 measurement points (determined by [1,max(2,round(n/40)),max(3,round(n/10)),max(4,round(n/2.5)),n]), in the present study, only the last measurement point (the final generation) was used. selfpenalty: g/snap by M+1 matrix that contains the local realised fitness cost of disallowing selfing in each deme ("(Npm - floor("number of produced gametes"/2))./Npm"). Please note that in Figs. B1-B3 in "The effect of the recombination rate between adaptive loci on the capacity of a population to expand its range" the outcrossing costs shown in Panel C are selfpenalty/(selfpenalty - 1). z: g/snap by M+1 matrix that contains the local trait mean in each in each deme. nohap: g/snap by M+1 matrix that contains the local number of unique haplotypes in each deme. meanfitness: g/snap by M+1 matrix that contains the local average Malthusian fitness in each deme. If desired the simulations can be prolonged, by launching the file "RangeExpansionCont.m" in the same folder as the output data from "RangeExpansion.m". Please modify TL and snap in RangeExpansionCont.m to set the total number of generations and the number of generations between measurement points. RangeExpansionCont.m uses DiploidMainCont.m as the main function, which in turn uses the same subfunction as DiploidMain.m (i.e. DiploidSelection.m, DiploidMutation.m, DiploidFertilisation.m and DiploidMigration.m).