ClearAll["Global`*"]

xMigration[k_, gf_, l_, t_, parameters_] := (1 - m) x[k, gf, l, t - 1, parameters] + m x[3 - k, gf, l, t - 1, parameters] /. parameters;

xMating[k_, gf1_, l1_, gf2_, l2_, t_, parameters_] := If[gf1 == 1 || gf1 == 3, If[gf2 == 1 || gf2 == 2, x[k, gf1, l1, t - 1, parameters]*xMigration[k, gf2, l2, t, parameters], 0], x[k, gf1, l1, t - 1, parameters]*xMigration[k, gf2, l2, t, parameters]]/Sum[xMigration[k, 1, L, t, parameters] + xMigration[k, 2, L, t, parameters] + KroneckerDelta[Mod[gf1, 2], 0] xMigration[k, 3, L, t, parameters] + KroneckerDelta[Mod[gf1, 2], 0] xMigration[k, 4, L, t, parameters], {L, 2}];


xSelection[k_, gf1_, l1_, gf2_, l2_, t_, parameters_] :=
  ((1 - (KroneckerDelta[l1, 3 - k] + KroneckerDelta[l2, 3 - k]) s) xMating[k, gf1, l1, gf2, l2, t, parameters])/Sum[Sum[Sum[Sum[(1 - (KroneckerDelta[L1, 3 - k] + KroneckerDelta[L2, 3 - k]) s) xMating[k, GF1, L1, GF2, L2, t, parameters], {GF1, 1, 4}], {L1, 1, 2}], {GF2, 1, 4}], {L2, 1, 2}] /. parameters

xRecombination[k_, gf1_, l1_, t_, parameters_] := Sum[If[GF1 == gf1 || GF2 == gf1,(*One of the chromosomes has to have the gametophytic factor haplotype gf1*)If[L1 == l1 || L2 == l1, (*One of the chromosomes has to have the trait allele l1*) If[GF1 == GF2 && L1 == L2, 1, (*If the chromosomes have the same genotype, all of their gametes will be this genotype*) If[GF1 == GF2 || L1 == L2, 1/2, (*If the chromosomes share an allele at the gf locus OR the trait locus, recombination doesn't matter and half of their gametes will be genotype gf l*)If[GF1 == gf1 && L1 == l1, 1/2 (1 - r),(*If one chromosome has the gamete genotype, half of nonrecombinant offspring will be be gf l*) If[GF2 == gf1 && L2 == l1, 1/2 (1 - r), (*If the other chromosome has the gamete genotype, half of nonrecombinant offspring will be gf l*)1/2 r (*If neither parents has genotype pi ti, it can only be made through recombination*)]]]], 0], 0] xSelection[k, GF1, L1, GF2, L2, t, parameters] /. parameters, {GF1, 4}, {L1, 2}, {GF2, 4}, {L2, 2}]

(*Memoization*)
x[k_, gf_, l_, t_, parameters_] := x[k, gf, l, t, parameters] = xRecombination[k, gf, l, t, parameters]

(*SL*)x[1, 1, 1, 0, parameters_] := 0.05;
(*ML*)x[1, 2, 1, 0, parameters_] := 0.95;
(*FL*)x[1, 3, 1, 0, parameters_] := 0;
(*gL*)x[1, 4, 1, 0, parameters_] := 0;
(*Sl*)x[1, 1, 2, 0, parameters_] := 0;
(*Ml*)x[1, 2, 2, 0, parameters_] := 0;
(*Fl*)x[1, 3, 2, 0, parameters_] := 0;
(*gl*)x[1, 4, 2, 0, parameters_] := 0;
(*SL*)x[2, 1, 1, 0, parameters_] := 0;
(*ML*)x[2, 2, 1, 0, parameters_] := 0;
(*FL*)x[2, 3, 1, 0, parameters_] := 0;
(*gL*)x[2, 4, 1, 0, parameters_] := 0;
(*Sl*)x[2, 1, 2, 0, parameters_] := 0;
(*Ml*)x[2, 2, 2, 0, parameters_] := 0;
(*Fl*)x[2, 3, 2, 0, parameters_] := 0;
(*gl*)x[2, 4, 2, 0, parameters_] := 1;

params = {m -> 0.01, s -> 0.5, r -> 0.01};

x[1, 1, 1, 10, params]//AbsoluteTiming