% MYERF Computation of erf(x) by Simpson's rule.

h=0.1;  x=0:h:3; 
Serf=zeros(1,31); 
for j=1:30
    mid=(x(j)+x(j+1))/2;
    temp=(h/6)*(exp(-x(j)^2) + 4*exp(-mid^2) + exp(-x(j+1)^2));
    Serf(j+1)=(2/sqrt(pi))*temp + Serf(j);
end
Serf
erf(x)
max(abs(Serf-erf(x)))