/* Integral(Midpoint Method) from a to b */
new; cls;
print midpointint(-100,0,1000000,&f);




proc midpointint(a,b,n,&f);
   local h,s,i,f:proc;
   h=(b-a)/n;
   s=0;
   i=1;
   do while i<=n;
      s=s+f(a+(2*i-1)/2*h);
      i=i+1;
   endo;
   s=h*s;
   retp(s);
endp;

fn f(x)=pdfn(x);