Hi, guys
I would like to know how to speed up numerical calculations of the inverse method in a game theoretical model Thank you very much!

Clear["`*"];
(*Parameter Definition*)a = 10;
c = {1, 2, 3, 8};
Kvals = {4, 4, 2, 1};
(*The fourth person's best response*)
q4Opt[q1_?NumericQ, q2_?NumericQ, q3_?NumericQ] := 
 Module[{q4}, 
  q4 /. FindMaximum[{(a - (q1 + q2 + q3 + q4)) q4 - c[[4]] q4, 
      0 <= q4 <= Kvals[[4]]}, {q4, 5}][[2]]]

(*The third person's best response*)
q3Opt[q1_?NumericQ, q2_?NumericQ] := 
 Module[{q3}, 
  q3 /. FindMaximum[{(a - (q1 + q2 + q3 + q4Opt[q1, q2, q3])) q3 - 
       c[[3]] q3, 0 <= q3 <= Kvals[[3]]}, {q3, 10}][[2]]]

(*The second person's best response*)
q2Opt[q1_?NumericQ] := 
 Module[{q2}, 
  q2 /. FindMaximum[{(a - (q1 + q2 + q3Opt[q1, q2] + 
            q4Opt[q1, q2, q3Opt[q1, q2]])) q2 - c[[2]] q2, 
      0 <= q2 <= Kvals[[2]]}, {q2, 12.5}][[2]]]

(*The first person's best decision*)
q1Opt = Module[{q1}, 
   q1 /. FindMaximum[{(a - (q1 + q2Opt[q1] + q3Opt[q1, q2Opt[q1]] + 
             q4Opt[q1, q2Opt[q1], q3Opt[q1, q2Opt[q1]]])) q1 - 
        c[[1]] q1, 0 <= q1 <= Kvals[[1]]}, {q1, 15}][[2]]];

(*Output*)
Print["Optimal decisions:"];
Print["q1 = ", q1Opt];
q2Val = q2Opt[q1Opt];
Print["q2 = ", q2Val];
q3Val = q3Opt[q1Opt, q2Val];
Print["q3 = ", q3Val];
q4Val = q4Opt[q1Opt, q2Val, q3Val];
Print["q4 = ", q4Val];