I'd go with an explicit summation using the trapezoidal rule. That way you need not figure out how to impose a specific sampling in NIntegrate. Can be done with Sum as below though certainly other variants are possible.

trapezoidal[func_, var_, p1_, p2_, n_] := With[{dx = (p2 - p1)/n},
  Sum[((func /. var -> N[j]) + (func /. var -> N[j - dx]))*dx/2, {j, 
    p1 + dx, p2, dx}]]

Examples:

Map[trapezoidal[Exp[-x^2], x, 0, 1, #] &, {20, 200, 2000}]
In[328]:= NIntegrate[Exp[-x^2], {x, 0, 1}]

(* Out[328]= 0.746824132812 *)

N[Integrate[Exp[-x^2], {x, 0, 1}]]

(* Out[334]= 0.7468241328121 *)
(* Out[333]= {0.74667083694, 0.74682259998, 0.746824117484} *)

Compare to NIntegrate and to N[] of the exact integral.

In[328]:= NIntegrate[Exp[-x^2], {x, 0, 1}]

(* Out[328]= 0.746824132812 *)

N[Integrate[Exp[-x^2], {x, 0, 1}]]

(* Out[334]= 0.7468241328121 *)

Be sure to acknowledge receiving assistance from Wolfram Community.

While you're at it, also try this:

WolframAlpha["Integrate e^(-(x^2)) from 0 to 1 using the Trapezoid Rule"]

This calls WolframAlpha and shows you what code you might use to solve this problem. You can see that it's possible to change the number of points used in the Trapezoid Rule.

There are more examples of this on Wolfram|Alpha's webpage:

http://www.wolframalpha.com/examples/NumericalIntegration.html

Here's how you can use Wolfram|Alpha to look at the 5 interval trapezoid rule for this integral:

http://www.wolframalpha.com/input/?i=5+interval+trapezoidal+rule+integrate+e%5E%28-%28x%5E2%29%29+on+%5B0%2C1%5D