program stab_euler
    logical :: top
    
    alpha = 2.

    tol = 1E-5
    !tol = 1E-3
    a = 0.
    b = 10.

    !h = 2E-2
    h = 1E-3
    !h = 1E-4
    !h = 1.15E-1

    do kk = 1, 5
    
       errSum = 0.0


       x = a
       top = .false.
       inode = 0

       ! IC
       y = exp(-alpha*(x-1)**2)

       open(11, file='u_h' , status='UNKNOWN')
       open(12, file='u_h_ref' , status='UNKNOWN')
       write(11,'(6es12.4)') h, x, y, exp(-alpha*(x-1)**2)

10     continue
       ynew = y + h * f(x, y, alpha)

       ypp = ( f(x+h, ynew, alpha) - f(x, y, alpha)) /h
       hnew = sqrt(2*tol / abs(ypp) )
       inode = inode + 1

       err_local = ypp*h*h / 2

       !if (err_local > tol) then
       !   ! step is refused
       !   write(12,'(60es12.4)') h, x+h, ynew, exp(-alpha*(x-1)**2), err_local, tol, ypp, hnew
       !   h = hnew * 0.9
       !   goto 10
       !endif

       errSum = errSum +   ( ( ynew - exp(-alpha*(x+h-1)**2))**2 &
            + ( y - exp(-alpha*(x-1)**2) )**2 ) / 2 * h


       x = x + h
       y = ynew

       write(11,'(6es12.4)') h, x, y, exp(-alpha*(x-1)**2), ypp, hnew

       if(.not. top .and. x >= 1.) then
          top = .true.
          errTop =  abs(y - exp(-alpha*(x-1)**2) )
       endif


       ! comment the following line for a fixed time step
       h = hnew

       if( x < b) goto 10
       errFin =  abs(y - exp(-alpha*(x-1)**2) )

       write(*,'(a40, 5es12.4,i8)') 'h, tol, errors at t= 1 and t=10, nodes :', &
            h, tol, errTop, errFin, sqrt(errSum), inode

       close(11)
       close(12)

       h = h/ 10
       tol = tol / 10
    enddo
  end program


  function f(x, y, alpha)
    real :: f
    real :: x, y, alpha

    f = -2* alpha * (x-1) * y
    
  end function f