proc iml;
   start myQntl(q, x, p);         /* definition 5 from UNIVARIATE doc */
      y = colvec(x);
      call sort(y,1);
      n = nrow(y);                   /* assume nonmissing data */
      q = j(ncol(p),1);
      do i = 1 to ncol(p);
         j  = (n+1)*p;            /* position in ordered data */
         j1 = floor(j); j2 = ceil(j);/* indices into ordered data */
         print j j1 j2;
         if j1=j2 then               /* return a datum */
            q = y[j1];
         else                        /* interpolate between data */
            q = (y[j2]+y[j1])/2;
      end;
   finish;
   /* test it */
   x = ranuni(j(100,1));
   p = {0.1 0.25 0.5 0.75 0.9};
   call MyQntl(q, x, p);
   print q;