function []=Gene_Dynamics(W,a,MaxT)
%
% Gene_Dynamics(W, a0, MaxTime)
%    where W is a 3x1 vector of the fitness of types aa, aA and AA.
%    a0 is the initial ratio of allele a in the population
%    MaxTime is the duration of the dynamics.
%

if nargin<3
    MaxT=100;
    if nargin<2
        a=1e-3;
        if nargin<1
            W=[2 1 1];
        end
    end
end

[T P]=ode45( @(t,y)diff_equations(t,y,W), [0 MaxT], a);

subplot(2,1,1);
plot(T,P,'k','linewidth',1);
xlabel 'Time'
ylabel 'Proportion of a allele in population'

subplot(2,1,2);
h=plot(T,P.*P,'-r',T,2*P.*(1-P),'-m',T,(1-P).*(1-P),'b'); set(h,'LineWidth',1); set(h(2),'Color',[0.7 0 0.8]);
xlabel 'Time'
ylabel 'Proportion of genotypes in population'
legend(h,'aa','aA','AA');


end

function [da]=diff_equations(t,a,W)

da=a*(1-a)*(W(1)*a + W(2)*(1-2*a) - W(3)*(1-a));

end