f[\[Theta]1_, \[Theta]2_] := 
  Module[{eqns, soln, A, B, plot, \[Alpha]2, d\[Theta]1 = 1, 
    d\[Theta]2 = 0, m1 = 1, m2 = 1, l1 = 1, l2 = 1, 
    g = 1, \[Alpha] = 0},
   \[Alpha]2 = \[Alpha] B'[t];
   eqns = {m2 l2 B''[t] + m2 l1 A''[t] Cos[A[t] - B[t]] - 
       m2 l1 (A'[t]^2) Sin[A[t] - B[t]] + 
       m2 g Sin[B[t]] + \[Alpha]2 m2 l2 == 
      0, (m1 + m2) l1 A''[t] + m2 l2 B''[t] Cos[A[t] - B[t]] + 
       m2 l2 (B'[t]^2) Sin[A[t] - B[t]] + 
       g (m1 + m2) Sin[A[t]] + \[Alpha] A'[t] (m1 + m2) == 0, 
     A[0] == \[Theta]1, A'[0] == d\[Theta]1, B[0] == \[Theta]2, 
     B'[0] == d\[Theta]2};
   soln = 
    NDSolve[eqns, {A, B}, {t, 0, 5}, MaxSteps -> Infinity, 
     PrecisionGoal -> 10];
   {Mod[First[A[5] /. soln], 2 Pi]/(2 Pi), 
    Mod[First[B[5] /. soln], 2 Pi]/(2 Pi)}
   ];
color[{a_, b_}] := Module[{c = ColorData["Rainbow"][a]},
  {c[[1]], c[[2]], c[[3]]} b
  ];
\[Delta] = 0.1;
xys = Table[{x, y}, {y, 0, 2 Pi, \[Delta]}, {x, 0, 2 Pi, \[Delta]}];
img = Map[color[f @@ #] &, xys, {2}];
Image[img]