program shri

    real(8) :: Y1(2) = 0d0,Y2(2) = 0d0, Eleft = 0d0,Eright = 0d0,E(20) = 0d0,Esearch = 0d0
    real(8) :: step = 0d0,acc = 1.d-5,h = 0d0,hmin = 1.d-10,xn = 0d0,xnd = 0d0,first = 0d0
    real(8) :: second = 0d0,third = 0d0,coeff = 0d0,u0 = 0d0,psi(201) = 0d0,norma = 0d0,x0 = 0d0
    real(8) :: slandau = 0d0, eps = 0d0,elandau(20) = 0d0
    integer :: j = 1,amount=999
    common /block/ Esearch, U0,coeff
    external pot

    coeff = 4.205208417d0

    !------
    U0=1.0d0
    !------

    step = U0 / 1.d3
    slandau = 1.d0/2.d0*(-1.d0+dsqrt(1.d0+coeff*4.d0*U0))

    !подсчет из Ландау-----------

    do while(eps.le.slandau)
        elandau(j) = -(1.d0/(4.d0*coeff))*(-(1.d0+2.d0*(eps))+dsqrt(1.d0+coeff*4.d0*U0))**2
        eps = eps + 1.
        j=j+1
    end do
    write(*,*)
    j=1
    !------------------------

    do i=1,amount

        Eleft=-u0+step*(i-1)
        Eright=-u0+step *i

        xnd=0.d0

        Esearch=Eleft
        xn=-200.d0
        Y1(1)=dexp(dsqrt(-coeff*Esearch)*xn)
        Y1(2)=Y1(1)*dsqrt(-coeff*Esearch)
        h=0.001d0
        call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
        xn=200.d0
        h=-0.0001d0
        Y2(1)=dexp(-dsqrt(-coeff*Esearch)*xn)
        Y2(2)=-Y2(1)*dsqrt(-coeff*Esearch)
        call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
        first=y1(2)/y1(1)-y2(2)/y2(1)
        Esearch=Eright
        xn=-100.d0
        Y1(1)=dexp(dsqrt(-coeff*Esearch)*xn)
        Y1(2)=Y1(1)*dsqrt(-coeff*Esearch)
        h=0.001d0
        call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
        xn=200.d0
        h=-0.0001d0
        Y2(1)=dexp(-dsqrt(-coeff*Esearch)*xn)
        Y2(2)=-Y2(1)*dsqrt(-coeff*Esearch)
        call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
        second=y1(2)/y1(1)-y2(2)/y2(1)
        if((first*second)<0) then
            do while(dabs(Eleft-Eright)>1.d-5)
                Esearch=(Eleft+Eright)/2
                E(j)=Esearch
                xn=-200.d0
                Y1(1)=dexp(dsqrt(-coeff*Esearch)*xn)
                Y1(2)=Y1(1)*dsqrt(-coeff*Esearch)
                h=0.1d-5
                call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
                xn=200.d0
                h=-0.1d-5
                Y2(1)=dexp(-dsqrt(-coeff*Esearch)*xn)
                Y2(2)=-Y2(1)*dsqrt(-coeff*Esearch)
                call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
                third=y1(2)/y1(1)-y2(2)/y2(1)
                if((first*third)<0) then
                    Eleft=Esearch
                else
                    Eright=Esearch
                end if
            end do
            j=j+1
        else
            continue
        end if
    end do

    write(*,*) 'Energy'
    write(*,*) 'Prog                    Theor'
    do i=1,j-1
        write(*,*) E(i), elandau(i), E(i)-elandau(i)
    end do
    write(*,*) 'levels: ', j-1

    open(12, file = 'd1.txt')
    Esearch=E(1)
    x0=-100.d0
    Y1(1)=dexp(dsqrt(-coeff*Esearch)*x0)
    Y1(2)=Y1(1)*dsqrt(-coeff*Esearch)
    step =0.1d0
    h=0.0001d0
    xn=-100.d0
    xnd=-10.d0
    call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
    x0=-10.d0
    acc=1.d-7
    do i=1,100
        Y2(1)=Y1(1)
        Y2(2)=Y1(2)
        xnd=-10.d0+step * i
        xn=x0
        call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
        psi(i)=Y2(1)
    end do
    x0=100.d0
    Y1(1)=dexp(-dsqrt(-coeff*Esearch)*x0)
    Y1(2)=-Y1(1)*dsqrt(-coeff*Esearch)
    s=0.1d0
    h=-0.0001d0
    xn=100.d0
    xnd=10.d0
    call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
    x0=10.d0
    acc=1.d-7
    do i=1,100
        Y2(1)=Y1(1)
        Y2(2)=Y1(2)
        xnd=10.d0-step *i
        xn=x0
        call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
        psi(201-i)=Y2(1)
    end do

    !norma=0.d0
    !do i=1,67
    !    norma=norma+step /3*((psi(3*(i-1)+1))**2+4*(psi(3*(i-1)+2))**2+(psi(3*(i-1)+3))**2)
    !end do
    !norma=dsqrt(norma)
    do i=1,201
        write(12,*) -10.d0+step*i, psi(i) !psi(i)/norma
    end do

    !--------------------------------------

    !open(13, file = 'd2.txt')
    !Esearch=E(2)
    !x0=-100.d0
    !Y1(1)=dexp(dsqrt(-coeff*Esearch)*x0)
    !Y1(2)=Y1(1)*dsqrt(-coeff*Esearch)
    !step =0.1d0
    !h=0.0001d0
    !xn=-100.d0
    !xnd=-10.d0
    !call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
    !x0=-10.d0
    !acc=1.d-7
    !do i=1,100
    !    Y2(1)=Y1(1)
    !    Y2(2)=Y1(2)
    !    xnd=-10.d0+step *i
    !    xn=x0
    !    call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
    !    psi(i)=Y2(1)
    !end do
    !x0=100.d0
    !Y1(1)=dexp(-dsqrt(-coeff*Esearch)*x0)
    !Y1(2)=-Y1(1)*dsqrt(-coeff*Esearch)
    !step =0.1d0
    !h=-0.0001d0
    !xn=100.d0
    !xnd=10.d0
    !call merson(xn,xnd,acc,h,hmin,2,Y1,pot)
    !x0=10.d0
    !acc=1.d-7
    !do i=1,100
    !    Y2(1)=Y1(1)
    !    Y2(2)=Y1(2)
    !    xnd=10.d0-s*i
    !    xn=x0
    !    call merson(xn,xnd,acc,h,hmin,2,Y2,pot)
    !    psi(201-i)=Y2(1)
    !end do

    !norma=0.d0
    !do i=1,67
    !    norma=norma+step /3*((psi(3*(i-1)+1))**2+4*(psi(3*(i-1)+2))**2+(psi(3*(i-1)+3))**2)
    !end do
    !norma=dsqrt(norma)
    !do i=1,201
    !    write(13,*) -10.d0+step *i, psi(i)/norma
    !end do

end

!для Мерсона в потенциале ch**2---------
subroutine pot(x, y, dy)
	real(8) :: x, y(2), dy(2), U0, Esearch, coeff
	common /block/ Esearch, U0,coeff
	dy(1)=y(2)
	dy(2)=-coeff*(Esearch+U0/(dcosh(x))**2)*y(1)
end