Can be used ParametricNDSolveValue

F = 1.2;
sol = ParametricNDSolveValue[{f'[
       x] (1 + 20 f[x] (1 - f[x]) f[x]^2) == -f[x], f[0] == a}, 
   f, {x, 0, 5}, {a}];
M[a_?NumericQ] := (NIntegrate[sol[a][x], {x, 0, 5}] - F)^2


sol1 = NMinimize[{M[a], 0 < a < 1}, a]
(*{1.03015*10^-21, {a -> 0.69524}}*)

Can be used ParametricNDSolveValue

F = 1.2;
sol = ParametricNDSolveValue[{f'[
       x] (1 + 20 f[x] (1 - f[x]) f[x]^2) == -f[x], f[0] == a}, 
   f, {x, 0, 5}, {a}];
M[a_?NumericQ] := (NIntegrate[sol[a][x], {x, 0, 5}] - F)^2


sol1 = NMinimize[{M[a], 0 < a < 1}, a]
(*{1.03015*10^-21, {a -> 0.69524}}*)

General view of the solution and optimal solution

{Plot[Evaluate[Table[sol[a][x], {a, 0.1, 1, .1}]], {x, 0, 5}], 
 Plot[sol[a][x] /. Last[sol1], {x, 0, 5}]}