      subroutine netstuf(x,w,nodes,ann,anx,anxx,gann,ganx,ganxx)
* Sigmoid Neural Network, with one hidden layer and linear output.
*    INPUT:
* x is the feature in ode space dimension
* w the vector of the weights (a total of 3*NODES)
* NODES the number of nodes
*    OUTPUT:
* ANN,ANX,ANXX: Values for the network and its first two x-derivatives
* GANN, GANX, GANXX: The gradients wrt w, of ANN, ANX, ANXX
      implicit double precision (a-h,o-z)
      dimension w(3*nodes),gann(3*nodes),ganx(3*nodes),ganxx(3*nodes)
      sig(z)=1./(1.+exp(-z))
      ann  = 0.d0
      anx  = 0.d0
      anxx = 0.d0
      do i=1,nodes
       ia = 3*i-2
       iw = 3*i-1
       ib = 3*i
       arg = w(iw)*x+w(ib)
       si = sig(arg)
       sp = si*(1.-si)
       sd = sp*(1.-2*si)
       st = sp*(1.-6*sp)
       ann = ann + w(ia)*si
       gann(ia)  = si
       gann(iw)  = w(ia)*x*sp
       gann(ib)  = w(ia)*sp
       anx = anx + w(ia)*w(iw)*sp
       ganx(ia)  = w(iw)*sp
       ganx(iw)  = w(ia)*(sp+w(iw)*x*sd)
       ganx(ib)  = w(ia)*w(iw)*sd
       anxx= anxx + w(ia)*w(iw)**2*sd
       ganxx(ia) = w(iw)**2*sd
       ganxx(iw) = w(ia)*(2*w(iw)*sd+w(iw)**2*x*st)
       ganxx(ib) = w(ia)*w(iw)**2*st
      enddo
      end
*
      function funmin(w,n)
* This is the relevant Loss function for solving the ODE:
!      y(x)'' + x**2*y(x) =(1+x**2)*exp(x), x in [0,1]
!      with y'(0)= 1  and y'(1)= e (e = exp(1.))
! using the trial solution: 
!   Y(x,w)= -0.5*(1+N(0,w))*(x-1)**2+0.5*(e-N(1,w))*x**2 +x*(x-1)*N(x,w)
*
      implicit double precision (a-h,o-z)
      save
      include 'PARA'
*     parameter(x0=0.d0,x1=1.d0, nop=50,nx=200)
      dimension w(n)
      dimension gann(nx),ganx(nx),ganxx(nx)
      common/tograd/ xg(nop)
      data itric/0/
      if (itric.eq.0) then
       itric = 1
       e=exp(1.d0)
       nodes = n/3
       if (3*nodes .ne. n) stop'N should be multiple of 3'
       if (n .gt. nx) stop 'increase NX >=N'
       do i=1,nop
         xg(i) = x0+(i-1.d0)*(x1-x0)/(nop-1.d0)
       enddo
      endif

      funmin = 0.d0
        x = 0.d0
        call netstuf(x,w,nodes,ann0,anx0,anxx0,gann,ganx,ganxx)
        x = 1.d0
        call netstuf(x,w,nodes,ann1,anx1,anxx1,gann,ganx,ganxx)
      do i = 1,nop
        x = xg(i)
        call netstuf(x,w,nodes,ann,anx,anxx,gann,ganx,ganxx)
        part1 = 2*ann+2*(2*x-1)*anx+x*(x-1)*anxx+e-ann1-1-ann0
        part2 =  x**2*(-0.5*(x-1)**2*(1+ann0)+0.5*(e-ann1)*x**2
     >          +x*(x-1)*ann)
        part3 = exp(x)*(1.+x**2)
        funmin = funmin + (part1+part2-part3)**2
      enddo
        funmin = funmin/nop
        call getfla(1,pri)
*  PRINTING for GNUPLOT
        if (nint(pri) .eq. 1) then
            iu = nunit()
            open(iu,file='COMPARE')
            call setfla(1,0.d0)
            x = 0.d0
            call netstuf(x,w,nodes,ann0,anx0,anxx0,gann,ganx,ganxx)
            x = 1.d0
            call netstuf(x,w,nodes,ann1,anx1,anxx1,gann,ganx,ganxx)
          do i = 1,nop
            x = xg(i)
            call netstuf(x,w,nodes,ann,anx,anxx,gann,ganx,ganxx)
            y = -0.5*(1+ann0)*(x-1)**2+0.5*(e-ann1)*x**2 +x*(x-1)*ann
            write(iu,*)x,y,exp(x),y-exp(x)
          enddo
          close(iu)
        endif
*  END of PRINTING
      end
      subroutine granal(n,w,grad)
      implicit double precision (a-h,o-z)
      include 'PARA'
      dimension w(n),grad(n)
      dimension gann(nx),ganx(nx),ganxx(nx),gann0(nx),gann1(nx)
      common/tograd/ xg(nop)
      nodes=n/3
      e=exp(1.d0)
      do j=1,n
        grad(j) =0.d0
      enddo
        x = 0.d0
        call netstuf(x,w,nodes,ann0,anx0,anxx0,gann0,ganx,ganxx)
        x = 1.d0
        call netstuf(x,w,nodes,ann1,anx1,anxx1,gann1,ganx,ganxx)
        do i = 1,nop
          x = xg(i)
          call netstuf(x,w,nodes,ann,anx,anxx,gann,ganx,ganxx)
          part1 = 2*ann+2*(2*x-1)*anx+x*(x-1)*anxx+e-ann1-1-ann0
          part2 =  x**2*(-0.5*(x-1)**2*(1+ann0)+0.5*(e-ann1)*x**2
     >            +x*(x-1)*ann)
          part3 = exp(x)*(1.+x**2)
          parta = part1+part2-part3
          do j = 1,n
            gart1 = 2*gann(j)+2*(2*x-1)*ganx(j)+x*(x-1)*ganxx(j)
     >             -(gann0(j)+gann1(j)) 
            gart2 = x**2*(-0.5*(x-1)**2*gann0(j)-0.5*x**2*gann1(j)
     >            + x*(x-1)*gann(j))
            grad(j) = grad(j) + 2*parta*(gart1+gart2)
          enddo
        enddo
      do j=1,n
        grad(j) = grad(j)/nop
      enddo

        end  






