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




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

fn f(x)=pdfn(x);