Решение диф. уравнения по формуле Симпсона

Simpson38[a0_, b0_, m0_] := 
  Module[{a = N[a0], b = N[b0], h, m = m0, X}, 
       h = (b - a)/(3*m); Subscript[X, k_] = a + k*h; 
        Return[((3*h)/8)*Sum[f[Subscript[X, 3*k - 3]] + 
                3*f[Subscript[X, 3*k - 2]] + 
       3*f[Subscript[X, 3*k - 1]] + 
                f[Subscript[X, 3*k]], {k, 1, m}]]; ]; 
f[x_, y_] = y*y*Exp[x] - 2*y; 
a = 0.; 
b = 5; 
M = {1, 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 
   35, 
       40}; 
tbl = Table[j, {j, 1, 20}]; 
Print[""]; 
For[j = 1, j <= 20, j++, 
     tbl[[j]] = {M[[j]], 
     NumberForm[Simpson38[a, b, M[[j]]], {8, 6}]}; ]; 
Print[TableForm[tbl, TableHeadings -> {None, {"m sample points", 
             "Simpson38[a,b,m]"}}, TableSpacing -> {2, 1}]]; 
Needs["Graphics`Colors`"]; 
Plot[f[x, y], {x, 0, 2}, PlotRange -> {{0, 2.01}, {0, 2.01}}, 
     Ticks -> {Range[0, 2, 0.5], Range[0, 2, 0.5]}, Axes -> True, 
     PlotStyle -> Magenta];