This is a code I developed a while ago. It uses Simpson method.
#include <iostream>
#include <math.h>
using namespace std;
double NormSDist(double x);
template <typename T> T SimpsonMethod(T (*pfunc)(T), float a, float b,int n=100);
int main()
{
double area=SimpsonMethod<double>(&NormSDist,-5,1.2);
cout << area << endl;
return 0;
}
template <typename T> T SimpsonMethod(T (*pfunc)(T), float a, float b,int n)
{
T h,x,y,retVal;
h=(b-a)/n;
x=a;
y=(*pfunc)(x);retVal=y;
for(int i=1;i<n;i++)
{
x=a+i*h;
if(i%2==0) y=2*(*pfunc)(x);
if(i%2==1) y=4*(*pfunc)(x);
retVal=retVal+y;
}
x=b;y=(*pfunc)(x);retVal=retVal+y;
return retVal*h/3;
}
double NormSDist(double x)
{
return 1/sqrt(2*3.1415926536)*exp(-0.5*x*x);
}
Hope it helps.