Program SIMPSON_2VAR
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DOUBLE PRECISION X2(1000000)

        a=0d0
        b=1d0
        N=10000
        h=(b-a)/N


        do i=1,N
        ! ahora llamamos a las abcissas finales x2
        x2(i)=a+i*h
        enddo

        sum1=0d0
        sum2=0d0
        nimp= n/2-1
        np = n/2

        ! El cambio aqui es la utilizacion de la funcion g2
        do i=1,nimp
        sum1=sum1+2*g2(x2(2*i))
        enddo

        do i=1,np
        sum2=sum2+4*g2(x2(2*i-1))
        enddo

        ! 3. Calculamos el resultado final en el programa principal
        
        r=h/3d0*(g2(a)+sum1+sum2+g2(b))
        
        open(20,File="ejerc10simp_comp_2var.xls")



        write(20,*)r
        call cpu_time(finish)
        write(20,30)finish-start
30      Format(e18.8)
        close(20)
        end
        
       ! 1. definimos la funcion de dos variables a integrar
        FUNCTION  g1(x1,x2)
        Implicit Double Precision (A-H,O-Z)
 !       parameter (a,b)

        g1=exp(-(x1*x1)-(x2*x2))*cos(2d0*x1*x2)

        END

       ! 2. Hacemos el primer paso de integracion reduciendo una
       ! variable , de modo analogo al programa simpson compuesto
       ! para una variable
       
        Function g2(x2)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DOUBLE PRECISION x1(1000000)
         N=10000
         sum1=0d0
         sum2=0d0
         nimp= n/2-1
         np = n/2
         a=0d0
         b=1d0
        h=(b-a)/N


        do i=1,N
        x1(i)=a+i*h
        enddo


        do i=1,nimp
        sum1=sum1+2*g1(x1(2*i),x2)

        enddo

        do i=1,np
        sum2=sum2+4*g1(x1(2*i-1),x2)
        enddo

       ! asignacion de g2
       
        g2=h/3d0*(g1(a,x2)+sum1+sum2+g1(b,x2))
        
        END