I am trying to fit the model to our flu outbreak data, by finding the parameter values $\beta$ and $\gamma$ that maximize the likelihood.

(*total size of population*) n = 763;

(*days*)
tmax = 14;

(*Flu dataset*)
numberinfected = {3, 8, 26, 76, 225, 298, 258, 233, 189, 128, 68, 29, 
   14, 4};

pointsPlot = 
 ListPlot[numberinfected, PlotStyle -> Purple, 
  PlotTheme -> "Detailed", 
  FrameLabel -> {"Time (days)", "Number Infected"}, 
  PlotLegends -> {"Total cases"}]

(*Model equations*)
ClearAll[\[Beta], \[Gamma]]
sireqns = {s'[t] == -\[Beta]*s[t]*i[t]/n, 
   i'[t] == \[Beta]*s[t]*i[t]/n - \[Gamma]*i[t], 
   r'[t] == \[Gamma]*i[t]};
(*Initial conditions*)
initialConditions = {s[0] == n-1, i[0] == 1, r[0] == 0};

(* Solving the differential equations *)
sol = ParametricNDSolve[{sireqns, initialConditions}, {s, i, r}, {t, 
    0, tmax}, {\[Beta], \[Gamma]}];

predicted = i /. sol;

(*Likelihood function*)
likelihood[\[Beta]_?NumericQ, \[Gamma]_?NumericQ] := 
  Exp[Total[(numberinfected - predicted)^2]];


(*Maximize the likelihood*)
result = 
  NMaximize[{Log[likelihood[\[Beta], \[Gamma]]], 0 <= \[Beta] <= 1, 
    0 <= \[Gamma] <= 1}, {\[Beta], \[Gamma]}];

optimalParameters = {\[Beta], \[Gamma]} /. Last[result]

I am looking for this solution

I tried to write some code similar to Writing a sum-of-squares function or SIR model fits but towards the end, I can't find my error. Any suggestions?