a = 0.736;
a0 = 0.5;
sol =
NDSolve[
{AA'[t] == a^2,
AB'[t] == a*(1-a),
BB'[t] == (1-a)^2,
AA[0] == a0,
BB[0] == 1-a0,
AB[0] == 0
},
{AA,AB,BB},
{t,40}
]

Plot[
{AA[t]/(AA[t]+AB[t]+BB[t])/.sol,
AB[t]/(AA[t]+AB[t]+BB[t])/.sol,
BB[t]/(AA[t]+AB[t]+BB[t])/.sol
},
{t,0,40},
PlotStyle->{Red,Purple,Blue}
]

f[a_] := a^2/(1-a(1-a));
f[0.736]
f[0.126]

a = 0.736;
p = 0;
a0 = 0.5;
sol =
NDSolve[
{AA'[t] == a*(1-p),
AB'[t] == a*p+ (1-a)*p,
BB'[t] == (1-a)*(1-p),
AA[0] == a0,
BB[0] == 1-a0,
AB[0] == 0
},
{AA,AB,BB},
{t,40}
]

Plot[
{AA[t]/(AA[t]+AB[t]+BB[t])/.sol,
AB[t]/(AA[t]+AB[t]+BB[t])/.sol,
BB[t]/(AA[t]+AB[t]+BB[t])/.sol
},
{t,0,40},
PlotStyle->{Red,Purple,Blue}
]

Clear[p]; Clear[q]; Clear[AA]; Clear[AB]; Clear[BB];
a = 0.736;
a0 = 0.5;
e0 = 0;
sol =
NDSolve[
{AA'[t] == a*(1-p[t]),
AB'[t] == a*p[t] + (1-a)*q[t],
BB'[t] == (1-a)*(1-q[t]),
p'[t] == (1-a-p[t])p[t],
q'[t] == (a-q[t])q[t],
AA[0] == a0-e0/2,
BB[0] == 1-a0-e0/2,
AB[0] == e0,
p[0] == 0.01,
q[0] == 0.01
},
{AA,AB,BB,p,q},
{t,40}
]

Plot[
{AA[t]/(AA[t]+AB[t]+BB[t])/.sol,
AB[t]/(AA[t]+AB[t]+BB[t])/.sol,
BB[t]/(AA[t]+AB[t]+BB[t])/.sol
},
{t,0,40},
PlotStyle->{Red,Purple,Blue}
]