How to quad("sin(a*sin(x))*sin(x)",0,pi/2)?

[Francesco Potorti`]

   > for (i in 1:4)
     {
       quadr ("sin(p1*sin(x)).*sin(x)", 0, pi/2,1.e-6,0,i).Q
     }
       0.691  
       0.906  
       0.533  
      -0.104  

   The above is a piece of Rlab code that does what you want (I
   think). The point is: I translated quadr from Matlab source. It is
   written by Kirill K. Pankratov (address@hidden)

I found the following in the matlab archives, form the
contribution/integration directory or some like that, if I am not
wrong.  They should be public domain.  I also have the sources.

===================================================================
Index for integration
===================================================================

nit: <directory>

===================================================================

 Numerical Integration Toolbox
 
 MATLAB Toolbox for 1-D, 2-D, and n-D Numerical Integration
 
 The original 1-D routines were obtained from NETLIB and were 
 written by
          Howard Wilson
          Department of Engineering Mechanics
          University of Alabama
          Box 870278
          Tuscaloosa, Alabama 35487-0278
          Phone 205 348-1617
          Email address: HWILSON @ UA1VM.UA.EDU
 
 The rest of the routines were written by
          Bryce Gardner
          Ray W. Herrick Laboratories
          Purdue University
          West Lafayette, IN 47906
          Phone: 317-494-0231
          Fax:  317-494-0787
          Email:  address@hidden
 
 These are the general purpose integration routines: 
 
        quadg.m      -- High accuracy replacement for QUAD and QUAD8 (1-D)
        quad2dg.m    -- 2-D integration over a rectangular region
        quad2dggen.m -- 2-D integration over a general region
        quadndg.m    -- n-D integration over a n-D hyper-rectangular region

===================================================================

misc: <directory>  Integration of math expressions

===================================================================

dquad.m:

  DQUAD     Numerically integrates an expression using QUAD.
            I = DQUAD('1/(1+x^3)',a,b) approximates,
 
                                      a
                                     /
                                    |     1    
                                    | ----------- dX
                                    |        3
                                    |  (1 + x  )
                                   /
                                   b

      D. Thomas  2/93
      address@hidden

dquad8.m:

  DQUAD     Numerically integrates an expression using QUAD8.
            I = DQUAD8('exp(-x^2)/(1+x^(1/2))',a,b) approximates,
 
                                   a
                                  /        2
                                 |  exp(- x )
                                 |  --------- dx
                                 |        1/2
                                /    1 + x
                                b
      D. Thomas  2/93
      address@hidden

trapz.m:

  a = trapz(y,x)

  find the area under the curve y = f(x) using the trapezoidal rule
   x = vector containing the independent variable data points
       the elements of x need not be equally spaced
   y = f(x) = vector of the function values corresponding to the
       elements in x

  a = the area under f(x) in the interval defined by x

 
 Written by: Duane Hanselman, University of Maine, (207)-581-2246