subroutine dosim(result,f,h,n)

      implicit double precision(a-h,o-z)

      common/lager/trapez
      common/delta/ bek,del
      common/res/res(650)

      real*8 f(n),h(n),fpp(650),bek(650),res1(650),res2(650)
      data fpp,bek/1300*0.0d0/

      sum=0.
      nm=n-1
      sumc=0.
      res1(1)=0.0
      res2(1)=0.0
      res2(2)=0.0
      do 1 j=2,n
      sum=sum+h(j)*(f(j)+f(j-1))
      res1(j)=h(j)*(f(j)+f(j-1))/2.
1     continue
      fpp(2)=2.*((f(3)-f(2))/h(3) - (f(2)-f(1))/h(2))/(h(2)+h(3))
      trapez=sum/2. 
      do 2 j=3,nm
      fpp(j)=2.*((f(j+1)-f(j))/h(j+1)-(f(j)-f(j-1))/h(j))/(h(j)+
     1  h(j+1))
        dp=h(j)+h(j+1)
        dm=h(j)+h(j-1)
        res2(j)=(dm*fpp(j)+dp*fpp(j-1))*h(j)**3/(dm+dp)
        sumc=sumc+(dm*fpp(j)+dp*fpp(j-1))*h(j)**3/(dm+dp)
2     continue

      endcor1=(fpp(2)+(h(2)+h(3))*(fpp(2)-fpp(3))/ 
     1  (h(2)+2.*h(3)+h(4)))*h(2)**3
      endcor=endcor1+(fpp(n-1)+(h(n)+h(n-1))*(fpp(n-1)-fpp(n-2))/
     1  (h(n)+2.*h(n-1)+h(n-2)))*h(n)**3
      sumc=sumc+endcor
      result=trapez-sumc/12.

      do j=2,n
        res(j)=res(j-1)+res1(j)-res2(j)/12.
      enddo

      return
      end 

************************************************************************