clear all

% This Matlab script is used to compare the performance of CRISPR-based 
% split gene drive systems with different component parts (i.e. Cas9 
% promoters and target genes).

% Define line styles and legend entries to be used during plotting
LnSt = {'-';'--';':';'-.'};
Lgnd = {'kmo^{sgRNAs} (sds3G1)','kmo^{sgRNAs} (bgcn)','w^{U6b-GDe} (nup50)','*kmo^{sgRNAs} (sds3G1)'};

for Case = 1:4 % Loop through the four different combinations to be studied
    
    % Parameter sets for each promoter/target combination
    
    % Here relative fitness parameters (omeM or omeF) relate to a reduced
    % degree of survival to adulthood. Fecundity/fertility reduction 
    % parameters (f) are used only for nup50-Cas9.
    if(Case==1) % sds3-Cas9 / kmo
        phiM  = 0.71; phiF  = 0.87;
        % (1) bbaa,   (2) bbAa,     (3) bbAA,     (4) Bbaa,     (5) BbAa,     (6) BbAA,     (7) BBaa,     (8) BBAa,     (9) BBAA
        omeM1 = 1.00; omeM2 = 1.00; omeM3 = 0.81; omeM4 = 1.00; omeM5 = 1.00; omeM6 = 0.81; omeM7 = 0.45; omeM8 = 0.45; omeM9 = 1-(1-0.45)-(1-0.81);
        omeF1 = 1.00; omeF2 = 1.00; omeF3 = 0.81; omeF4 = 1.00; omeF5 = 1.00; omeF6 = 0.81; omeF7 = 0.91; omeF8 = 0.91; omeF9 = 1-(1-0.91)-(1-0.81);
        f1    = 1.00; f2    = 1.00; f3    = 1.00; f4    = 1.00; f5    = 1.00; f6    = 1.00; f7    = 1.00; f8    = 1.00; f9    = 1.00;
    elseif(Case==2) % bgcn-Cas9 / kmo
        phiM  = 0.36;  phiF = 0.54;
        omeM1 = 1.00; omeM2 = 1.00; omeM3 = 0.81; omeM4 = 1.00; omeM5 = 1.00; omeM6 = 0.81; omeM7 = 0.79; omeM8 = 0.79; omeM9 = 1-(1-0.79)-(1-0.81);
        omeF1 = 1.00; omeF2 = 1.00; omeF3 = 0.81; omeF4 = 1.00; omeF5 = 1.00; omeF6 = 0.81; omeF7 = 0.79; omeF8 = 0.79; omeF9 = 1-(1-0.79)-(1-0.81);
        f1    = 1.00; f2    = 1.00; f3    = 1.00; f4    = 1.00; f5    = 1.00; f6    = 1.00; f7    = 1.00; f8    = 1.00; f9    = 1.00;
    elseif(Case==3) % nup50-Cas9 / w
        phiM  = 0.37; phiF  = 0.63;
        omeM1 = 1.00; omeM2 = 1.00; omeM3 = 0.90; omeM4 = 1.00;  omeM5 = 1.00;  omeM6 = 0.90;  omeM7 = 1.00;  omeM8 = 1.00;  omeM9 = 0.90;
        omeF1 = 1.00; omeF2 = 1.00; omeF3 = 0.90; omeF4 = 1.00;  omeF5 = 1.00;  omeF6 = 0.90;  omeF7 = 1.00;  omeF8 = 1.00;  omeF9 = 0.90;
        f1    = 1.00; f2    = 1.00; f3    = 1.00; f4    = 0.922; f5    = 0.922; f6    = 0.922; f7    = 0.844; f8    = 0.844; f9    = 0.844;
    elseif(Case==4) % sds3-Cas9 / kmo (with w fitness)
        phiM  = 0.71; phiF  = 0.87;
        omeM1 = 1.00; omeM2 = 1.00; omeM3 = 0.90; omeM4 = 1.00; omeM5 = 1.00; omeM6 = 0.90; omeM7 = 0.45; omeM8 = 0.45; omeM9 = 1-(1-0.45)-(1-0.90);
        omeF1 = 1.00; omeF2 = 1.00; omeF3 = 0.90; omeF4 = 1.00; omeF5 = 1.00; omeF6 = 0.90; omeF7 = 0.91; omeF8 = 0.91; omeF9 = 1-(1-0.91)-(1-0.90);
        f1    = 1.00; f2    = 1.00; f3    = 1.00; f4    = 1.00; f5    = 1.00; f6    = 1.00; f7    = 1.00; f8    = 1.00; f9    = 1.00;
    end

    % Initial conditions (used for all four cases)
    % This represents a setup with 50% W-T males and 50% BbAa females
    % (1) bbaa,      (2) bbAa,      (3) bbAA,      (4) Bbaa,      (5) BbAa,      (6) BbAA,      (7) BBaa,      (8) BBAa,      (9) BBAA
      M(1,1) = 0.50; M(1,2) = 0.00; M(1,3) = 0.00; M(1,4) = 0.00; M(1,5) = 0.00; M(1,6) = 0.00; M(1,7) = 0.00; M(1,8) = 0.00; M(1,9) = 0.00;
      F(1,1) = 0.00; F(1,2) = 0.00; F(1,3) = 0.00; F(1,4) = 0.00; F(1,5) = 0.50; F(1,6) = 0.00; F(1,7) = 0.00; F(1,8) = 0.00; F(1,9) = 0.00;

    for i=2:51 % Loop through desired number of generations
        % Calculate proportional frequencies for each genotype and sex
        M1e = (omeM1/2)*(M(i-1,1)*F(i-1,1)*f1 + 0.5*M(i-1,1)*F(i-1,2)*f2 + 0.5*M(i-1,1)*F(i-1,4)*f4 + 0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,2)*F(i-1,1)*f1 + 0.25*M(i-1,2)*F(i-1,2)*f2 + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,1)*f1 + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.25*M(i-1,4)*F(i-1,4)*f4 + 0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF));
        M2e = (omeM2/2)*(0.5*M(i-1,1)*F(i-1,2)*f2 + M(i-1,1)*F(i-1,3)*f3 + 0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,1)*F(i-1,6)*f6 + 0.5*M(i-1,2)*F(i-1,1)*f1 + 0.5*M(i-1,2)*F(i-1,2)*f2 + 0.5*M(i-1,2)*F(i-1,3)*f3 + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.25*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + M(i-1,3)*F(i-1,1)*f1 + 0.5*M(i-1,3)*F(i-1,2)*f2 + 0.5*M(i-1,3)*F(i-1,4)*f4 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.5*M(i-1,4)*F(i-1,3)*f3 + 0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,4)*F(i-1,6)*f6 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.5*M(i-1,6)*F(i-1,1)*f1 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.25*M(i-1,6)*F(i-1,4)*f4 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF)+0.5*M(i-1,1)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,4)*F(i-1,5)*f5*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,5)*F(i-1,1)*f1*(phiM) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM)+ 0.25*M(i-1,5)*F(i-1,4)*f4*(phiM) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.125*(1-phiM)*M(i-1,5)*(phiF)*F(i-1,8)*f8);
        M3e = (omeM3/2)*(0.25*M(i-1,2)*F(i-1,2)*f2 + 0.5*M(i-1,2)*F(i-1,3)*f3 + 0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + 0.5*M(i-1,3)*F(i-1,2)*f2 + M(i-1,3)*F(i-1,3)*f3 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,3)*F(i-1,6)*f6 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.5*M(i-1,6)*F(i-1,3)*f3 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,6)*F(i-1,6)*f6+0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,3)*F(i-1,5)*f5*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.25*M(i-1,6)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,5)*F(i-1,3)*f3*(phiM)+ 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.25*M(i-1,5)*F(i-1,6)*f6*(phiM) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(phiF) + 0.125*(1-phiM)*M(i-1,5)*(phiF)*F(i-1,8)*f8);
        M4e = (omeM4/2)*(0.5*M(i-1,1)*F(i-1,4)*f4 + 0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + M(i-1,1)*F(i-1,7)*f7 + 0.5*M(i-1,1)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,2)*F(i-1,7)*f7 + 0.25*M(i-1,2)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,4)*F(i-1,1)*f1 + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.5*M(i-1,4)*F(i-1,4)*f4 + 0.25*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,7)*f7 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + M(i-1,7)*F(i-1,1)*f1 + 0.5*M(i-1,7)*F(i-1,2)*f2 + 0.5*M(i-1,7)*F(i-1,4)*f4 + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,1)*f1 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF));
        M5e = (omeM5/2)*(0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,1)*F(i-1,6)*f6 + 0.5*M(i-1,1)*F(i-1,8)*f8*(1-phiF) + M(i-1,1)*F(i-1,9)*f9 + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.25*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + 0.5*M(i-1,2)*F(i-1,7)*f7 + 0.5*M(i-1,2)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,2)*F(i-1,9)*f9 + 0.5*M(i-1,3)*F(i-1,4)*f4 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + M(i-1,3)*F(i-1,7)*f7 + 0.5*M(i-1,3)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.5*M(i-1,4)*F(i-1,3)*f3 + 0.25*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,6)*f6 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,4)*F(i-1,9)*f9 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.5*M(i-1,6)*F(i-1,1)*f1 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.5*M(i-1,6)*F(i-1,4)*f4 + 0.25*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,6)*F(i-1,7)*f7 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,7)*F(i-1,2)*f2 + M(i-1,7)*F(i-1,3)*f3 + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,7)*F(i-1,6)*f6 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,1)*f1 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,2)*f2 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,3)*f3 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) +0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + M(i-1,9)*F(i-1,1)*f1 + 0.5*M(i-1,9)*F(i-1,2)*f2 + 0.5*M(i-1,9)*F(i-1,4)*f4 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF)+0.5*M(i-1,1)*F(i-1,5)*f5*(phiF) + M(i-1,1)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,2)*F(i-1,8)*f8*(phiF) + 0.5*M(i-1,4)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,4)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM)+ 0.25*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.5*M(i-1,7)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,5)*F(i-1,1)*f1*(phiM) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,5)*F(i-1,4)*f4*(phiM)+ 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,7)*f7*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + M(i-1,8)*F(i-1,1)*f1*(phiM) + 0.5*M(i-1,8)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,8)*F(i-1,4)*f4*(phiM)+ 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF));
        M6e = (omeM6/2)*(0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + 0.25*M(i-1,2)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,2)*F(i-1,9)*f9 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,3)*F(i-1,6)*f6 + 0.5*M(i-1,3)*F(i-1,8)*f8*(1-phiF) + M(i-1,3)*F(i-1,9)*f9 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.5*M(i-1,6)*F(i-1,3)*f3 + 0.25*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,6)*F(i-1,6)*f6 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,6)*F(i-1,9)*f9 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,2)*f2 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,3)*f3 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + 0.5*M(i-1,9)*F(i-1,2)*f2 + M(i-1,9)*F(i-1,3)*f3 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,9)*F(i-1,6)*f6+0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,2)*F(i-1,8)*f8*(phiF) + 0.5*M(i-1,3)*F(i-1,5)*f5*(phiF) + M(i-1,3)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.5*M(i-1,6)*F(i-1,5)*f5*(phiF)+ 0.5*M(i-1,6)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,9)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,5)*F(i-1,3)*f3*(phiM) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF)+ 0.5*M(i-1,5)*F(i-1,6)*f6*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,9)*f9*(phiM) + 0.5*M(i-1,8)*F(i-1,2)*f2*(phiM) + M(i-1,8)*F(i-1,3)*f3*(phiM) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF)+ 0.5*M(i-1,8)*F(i-1,6)*f6*(phiM) + 0.5*M(i-1,5)*F(i-1,5)*f5*(phiM)*(phiF) + 0.5*M(i-1,5)*F(i-1,8)*f8*(phiM)*(phiF) + 0.5*M(i-1,8)*F(i-1,5)*f5*(phiM)*(phiF));
        M7e = (omeM7/2)*(0.25*M(i-1,4)*F(i-1,4)*f4 + 0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,7)*f7 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,7)*F(i-1,4)*f4 + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + M(i-1,7)*F(i-1,7)*f7 + 0.5*M(i-1,7)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,7)*f7 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,8)*f8*(1-phiF));
        M8e = (omeM8/2)*(0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,4)*F(i-1,6)*f6 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,4)*F(i-1,9)*f9 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.25*M(i-1,6)*F(i-1,4)*f4 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,6)*F(i-1,7)*f7 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,7)*F(i-1,6)*f6 + 0.5*M(i-1,7)*F(i-1,8)*f8*(1-phiF) + M(i-1,7)*F(i-1,9)*f9 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,7)*f7 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,9)*f9 + 0.5*M(i-1,9)*F(i-1,4)*f4 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF) + M(i-1,9)*F(i-1,7)*f7 + 0.5*M(i-1,9)*F(i-1,8)*f8*(1-phiF)+0.25*M(i-1,4)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,4)*F(i-1,8)*f8*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.125*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.5*M(i-1,7)*F(i-1,5)*f5*(phiF) + M(i-1,7)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM)+ 0.5*M(i-1,8)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.25*M(i-1,5)*F(i-1,4)*f4*(phiM) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,7)*f7*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + 0.5*M(i-1,8)*F(i-1,4)*f4*(phiM)+ 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF) + M(i-1,8)*F(i-1,7)*f7*(phiM) + 0.5*M(i-1,8)*F(i-1,8)*f8*(phiM)*(1-phiF));
        M9e = (omeM9/2)*(0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,6)*F(i-1,6)*f6 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,6)*F(i-1,9)*f9 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,9)*f9 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,9)*F(i-1,6)*f6 + 0.5*M(i-1,9)*F(i-1,8)*f8*(1-phiF) + M(i-1,9)*F(i-1,9)*f9+0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.125*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.25*M(i-1,6)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,6)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,8)*F(i-1,8)*f8*(phiF)*(1-phiM)+ 0.5*M(i-1,9)*F(i-1,5)*f5*(phiF) + M(i-1,9)*F(i-1,8)*f8*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.25*M(i-1,5)*F(i-1,6)*f6*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,9)*f9*(phiM)+ 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.5*M(i-1,8)*F(i-1,6)*f6*(phiM) + 0.5*M(i-1,8)*F(i-1,8)*f8*(phiM)*(1-phiF) + M(i-1,8)*F(i-1,9)*f9*(phiM) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(phiF) + 0.5*M(i-1,5)*F(i-1,8)*f8*(phiM)*(phiF)+ 0.5*M(i-1,8)*F(i-1,5)*f5*(phiM)*(phiF) + M(i-1,8)*F(i-1,8)*f8*(phiM)*(phiF));
        F1e = (omeF1/2)*(M(i-1,1)*F(i-1,1)*f1 + 0.5*M(i-1,1)*F(i-1,2)*f2 + 0.5*M(i-1,1)*F(i-1,4)*f4 + 0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,2)*F(i-1,1)*f1 + 0.25*M(i-1,2)*F(i-1,2)*f2 + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,1)*f1 + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.25*M(i-1,4)*F(i-1,4)*f4 + 0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF));
        F2e = (omeF2/2)*(0.5*M(i-1,1)*F(i-1,2)*f2 + M(i-1,1)*F(i-1,3)*f3 + 0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,1)*F(i-1,6)*f6 + 0.5*M(i-1,2)*F(i-1,1)*f1 + 0.5*M(i-1,2)*F(i-1,2)*f2 + 0.5*M(i-1,2)*F(i-1,3)*f3 + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.25*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + M(i-1,3)*F(i-1,1)*f1 + 0.5*M(i-1,3)*F(i-1,2)*f2 + 0.5*M(i-1,3)*F(i-1,4)*f4 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.5*M(i-1,4)*F(i-1,3)*f3 + 0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,4)*F(i-1,6)*f6 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.5*M(i-1,6)*F(i-1,1)*f1 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.25*M(i-1,6)*F(i-1,4)*f4 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF)+0.5*M(i-1,1)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,4)*F(i-1,5)*f5*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,5)*F(i-1,1)*f1*(phiM) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM)+ 0.25*M(i-1,5)*F(i-1,4)*f4*(phiM) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.125*(1-phiM)*M(i-1,5)*(phiF)*F(i-1,8)*f8);
        F3e = (omeF3/2)*(0.25*M(i-1,2)*F(i-1,2)*f2 + 0.5*M(i-1,2)*F(i-1,3)*f3 + 0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + 0.5*M(i-1,3)*F(i-1,2)*f2 + M(i-1,3)*F(i-1,3)*f3 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,3)*F(i-1,6)*f6 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.5*M(i-1,6)*F(i-1,3)*f3 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,6)*F(i-1,6)*f6+0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,3)*F(i-1,5)*f5*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.25*M(i-1,6)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,5)*F(i-1,3)*f3*(phiM)+ 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.25*M(i-1,5)*F(i-1,6)*f6*(phiM) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(phiF) + 0.125*(1-phiM)*M(i-1,5)*(phiF)*F(i-1,8)*f8);
        F4e = (omeF4/2)*(0.5*M(i-1,1)*F(i-1,4)*f4 + 0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + M(i-1,1)*F(i-1,7)*f7 + 0.5*M(i-1,1)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,2)*F(i-1,7)*f7 + 0.25*M(i-1,2)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,4)*F(i-1,1)*f1 + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.5*M(i-1,4)*F(i-1,4)*f4 + 0.25*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,7)*f7 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + M(i-1,7)*F(i-1,1)*f1 + 0.5*M(i-1,7)*F(i-1,2)*f2 + 0.5*M(i-1,7)*F(i-1,4)*f4 + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,1)*f1 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF));
        F5e = (omeF5/2)*(0.25*M(i-1,1)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,1)*F(i-1,6)*f6 + 0.5*M(i-1,1)*F(i-1,8)*f8*(1-phiF) + M(i-1,1)*F(i-1,9)*f9 + 0.25*M(i-1,2)*F(i-1,4)*f4 + 0.25*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + 0.5*M(i-1,2)*F(i-1,7)*f7 + 0.5*M(i-1,2)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,2)*F(i-1,9)*f9 + 0.5*M(i-1,3)*F(i-1,4)*f4 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + M(i-1,3)*F(i-1,7)*f7 + 0.5*M(i-1,3)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,4)*F(i-1,2)*f2 + 0.5*M(i-1,4)*F(i-1,3)*f3 + 0.25*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,6)*f6 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,4)*F(i-1,9)*f9 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,1)*f1 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.5*M(i-1,6)*F(i-1,1)*f1 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.5*M(i-1,6)*F(i-1,4)*f4 + 0.25*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,6)*F(i-1,7)*f7 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,7)*F(i-1,2)*f2 + M(i-1,7)*F(i-1,3)*f3 + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,7)*F(i-1,6)*f6 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,1)*f1 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,2)*f2 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,3)*f3 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) +0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + M(i-1,9)*F(i-1,1)*f1 + 0.5*M(i-1,9)*F(i-1,2)*f2 + 0.5*M(i-1,9)*F(i-1,4)*f4 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF)+0.5*M(i-1,1)*F(i-1,5)*f5*(phiF) + M(i-1,1)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,2)*F(i-1,8)*f8*(phiF) + 0.5*M(i-1,4)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,4)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM)+ 0.25*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.5*M(i-1,7)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,5)*F(i-1,1)*f1*(phiM) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,5)*F(i-1,4)*f4*(phiM)+ 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,7)*f7*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + M(i-1,8)*F(i-1,1)*f1*(phiM) + 0.5*M(i-1,8)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,8)*F(i-1,4)*f4*(phiM)+ 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF));
        F6e = (omeF6/2)*(0.125*M(i-1,2)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,2)*F(i-1,6)*f6 + 0.25*M(i-1,2)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,2)*F(i-1,9)*f9 + 0.25*M(i-1,3)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,3)*F(i-1,6)*f6 + 0.5*M(i-1,3)*F(i-1,8)*f8*(1-phiF) + M(i-1,3)*F(i-1,9)*f9 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,2)*f2 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,3)*f3 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.25*M(i-1,6)*F(i-1,2)*f2 + 0.5*M(i-1,6)*F(i-1,3)*f3 + 0.25*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,6)*F(i-1,6)*f6 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,6)*F(i-1,9)*f9 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,2)*f2 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,3)*f3 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + 0.5*M(i-1,9)*F(i-1,2)*f2 + M(i-1,9)*F(i-1,3)*f3 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,9)*F(i-1,6)*f6+0.25*M(i-1,2)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,2)*F(i-1,8)*f8*(phiF) + 0.5*M(i-1,3)*F(i-1,5)*f5*(phiF) + M(i-1,3)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.5*M(i-1,6)*F(i-1,5)*f5*(phiF)+ 0.5*M(i-1,6)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,9)*F(i-1,5)*f5*(phiF) + 0.25*M(i-1,5)*F(i-1,2)*f2*(phiM) + 0.5*M(i-1,5)*F(i-1,3)*f3*(phiM) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF)+ 0.5*M(i-1,5)*F(i-1,6)*f6*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,9)*f9*(phiM) + 0.5*M(i-1,8)*F(i-1,2)*f2*(phiM) + M(i-1,8)*F(i-1,3)*f3*(phiM) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF)+ 0.5*M(i-1,8)*F(i-1,6)*f6*(phiM) + 0.5*M(i-1,5)*F(i-1,5)*f5*(phiM)*(phiF) + 0.5*M(i-1,5)*F(i-1,8)*f8*(phiM)*(phiF) + 0.5*M(i-1,8)*F(i-1,5)*f5*(phiM)*(phiF));
        F7e = (omeF7/2)*(0.25*M(i-1,4)*F(i-1,4)*f4 + 0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,4)*F(i-1,7)*f7 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,7)*F(i-1,4)*f4 + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + M(i-1,7)*F(i-1,7)*f7 + 0.5*M(i-1,7)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,7)*f7 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,8)*f8*(1-phiF));
        F8e = (omeF8/2)*(0.125*M(i-1,4)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,4)*F(i-1,6)*f6 + 0.25*M(i-1,4)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,4)*F(i-1,9)*f9 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,4)*f4 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,7)*f7 + 0.25*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.25*M(i-1,6)*F(i-1,4)*f4 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,6)*F(i-1,7)*f7 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,7)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,7)*F(i-1,6)*f6 + 0.5*M(i-1,7)*F(i-1,8)*f8*(1-phiF) + M(i-1,7)*F(i-1,9)*f9 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,4)*f4 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,7)*f7 + 0.5*M(i-1,8)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,9)*f9 + 0.5*M(i-1,9)*F(i-1,4)*f4 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF) + M(i-1,9)*F(i-1,7)*f7 + 0.5*M(i-1,9)*F(i-1,8)*f8*(1-phiF)+0.25*M(i-1,4)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,4)*F(i-1,8)*f8*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.125*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.5*M(i-1,7)*F(i-1,5)*f5*(phiF) + M(i-1,7)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM)+ 0.5*M(i-1,8)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.25*M(i-1,5)*F(i-1,4)*f4*(phiM) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,7)*f7*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + 0.5*M(i-1,8)*F(i-1,4)*f4*(phiM)+ 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF) + M(i-1,8)*F(i-1,7)*f7*(phiM) + 0.5*M(i-1,8)*F(i-1,8)*f8*(phiM)*(1-phiF));
        F9e = (omeF9/2)*(0.0625*M(i-1,5)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.125*M(i-1,5)*(1-phiM)*F(i-1,6)*f6 + 0.125*M(i-1,5)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.25*M(i-1,5)*(1-phiM)*F(i-1,9)*f9 + 0.125*M(i-1,6)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,6)*F(i-1,6)*f6 + 0.25*M(i-1,6)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,6)*F(i-1,9)*f9 + 0.125*M(i-1,8)*(1-phiM)*F(i-1,5)*f5*(1-phiF) + 0.25*M(i-1,8)*(1-phiM)*F(i-1,6)*f6 + 0.25*M(i-1,8)*(1-phiM)*F(i-1,8)*f8*(1-phiF) + 0.5*M(i-1,8)*(1-phiM)*F(i-1,9)*f9 + 0.25*M(i-1,9)*F(i-1,5)*f5*(1-phiF) + 0.5*M(i-1,9)*F(i-1,6)*f6 + 0.5*M(i-1,9)*F(i-1,8)*f8*(1-phiF) + M(i-1,9)*F(i-1,9)*f9+0.125*M(i-1,5)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.125*M(i-1,5)*F(i-1,8)*f8*(phiF)*(1-phiM) + 0.25*M(i-1,6)*F(i-1,5)*f5*(phiF) + 0.5*M(i-1,6)*F(i-1,8)*f8*(phiF) + 0.25*M(i-1,8)*F(i-1,5)*f5*(phiF)*(1-phiM) + 0.5*M(i-1,8)*F(i-1,8)*f8*(phiF)*(1-phiM)+ 0.5*M(i-1,9)*F(i-1,5)*f5*(phiF) + M(i-1,9)*F(i-1,8)*f8*(phiF) + 0.125*M(i-1,5)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.25*M(i-1,5)*F(i-1,6)*f6*(phiM) + 0.25*M(i-1,5)*F(i-1,8)*f8*(phiM)*(1-phiF) + 0.5*M(i-1,5)*F(i-1,9)*f9*(phiM)+ 0.25*M(i-1,8)*F(i-1,5)*f5*(phiM)*(1-phiF) + 0.5*M(i-1,8)*F(i-1,6)*f6*(phiM) + 0.5*M(i-1,8)*F(i-1,8)*f8*(phiM)*(1-phiF) + M(i-1,8)*F(i-1,9)*f9*(phiM) + 0.25*M(i-1,5)*F(i-1,5)*f5*(phiM)*(phiF) + 0.5*M(i-1,5)*F(i-1,8)*f8*(phiM)*(phiF)+ 0.5*M(i-1,8)*F(i-1,5)*f5*(phiM)*(phiF) + M(i-1,8)*F(i-1,8)*f8*(phiM)*(phiF));
        % Calculate normalising factor (a.k.a. average fitness of population)
    	OmeBar = M1e+M2e+M3e+M4e+M5e+M6e+M7e+M8e+M9e+F1e+F2e+F3e+F4e+F5e+F6e+F7e+F8e+F9e;
        % Normalise results
    	M(i,1) = M1e/OmeBar; M(i,2) = M2e/OmeBar; M(i,3) = M3e/OmeBar; M(i,4) = M4e/OmeBar; M(i,5) = M5e/OmeBar; M(i,6) = M6e/OmeBar; M(i,7) = M7e/OmeBar; M(i,8) = M8e/OmeBar; M(i,9) = M9e/OmeBar;
        F(i,1) = F1e/OmeBar; F(i,2) = F2e/OmeBar; F(i,3) = F3e/OmeBar; F(i,4) = F4e/OmeBar; F(i,5) = F5e/OmeBar; F(i,6) = F6e/OmeBar; F(i,7) = F7e/OmeBar; F(i,8) = F8e/OmeBar; F(i,9) = F9e/OmeBar;
    end

    % Plotting Commands
    figure(1); hold on; box on;
    plot([0:i-1],M(:,4)+M(:,5)+M(:,6)+M(:,7)+M(:,8)+M(:,9)+F(:,4)+F(:,5)+F(:,6)+F(:,7)+F(:,8)+F(:,9),'r','LineStyle',LnSt{Case},'LineWidth',2,'DisplayName',Lgnd{Case}); % B
    plot([0:i-1],M(:,2)+M(:,3)+M(:,5)+M(:,6)+M(:,8)+M(:,9)+F(:,2)+F(:,3)+F(:,5)+F(:,6)+F(:,8)+F(:,9),'b','LineStyle',LnSt{Case},'LineWidth',2,'DisplayName',Lgnd{Case}); % A
end
% Add axis labels and select various appearance options
xlabel('Generation');
ylabel('Carrier Frequency');
legend('Location','best');
axis([0 i-1 0 1]);
set(gca,'FontName','Arial','FontSize',12);