Численное решение уравнения второго порядка

program solverdiffeq
integer, parameter:: n=100000
real dx,eps,a,b,tg0,yn,yk,dtg
integer i,k
real:: x(n),y(n),alfa(n),beta(n),y1(n),y2(n),f(n)
pi=3.1415926535897932384626433832795
a=0.0
b=2.*pi
tg0=0.0
dtg=0.1
yn=0.0
yk=0.0
eps=10.0**(-4)
eps1=1.0
dx=(b-a)/n
x(1)=a
do i=1,(n-1)
   x(i+1)=x(i)+dx
   alfa(i)=4*sin(x(i))
   beta(i)=cos(x(i))**(2.0)
end do
alfa(n)=4.*sin(x(n))
beta(n)=cos(x(n))**(2.0)
y(1)=yn
y(2)=y(1)+tg0*dx
k=0
do while (eps1>=eps)
         y(1)=yn
         y(2)=y(1)+(tg0-dtg)*dx
         y1(1)=yn
         y1(2)=y1(1)+tg0*dx
         y2(1)=yn
         y2(2)=y(1)+(tg0+dtg)*dx
         do i=2,(n-1)
             y(i+1)=beta(i)*dx**2+(2.-alfa(i)*dx**2)*y(i) - y(i-1)
            y1(i+1)=beta(i)*dx**2+(2.-alfa(i)*dx**2)*y1(i)-y1(i-1)
            y2(i+1)=beta(i)*dx**2+(2.-alfa(i)*dx**2)*y2(i)-y2(i-1)
         end do
        if ((abs(y(n)-yk)<abs(y1(n)-yk)) .and. (abs(y(n)-yk)<abs(y2(n)-yk))) then
            tg0=tg0-dtg
            eps1=abs(y(n)-yk)
            f=y
        end if
        if ((abs(y1(n)-yk)<=abs(y(n)-yk)) .and. (abs(y1(n)-yk)<=abs(y2(n)-yk))) then
            dtg=dtg/2.0
            eps1=abs(y1(n)-yk)
            f=y1
        end if
        if ((abs(y2(n)-yk)<abs(y(n)-yk)) .and. (abs(y2(n)-yk)<abs(y1(n)-yk))) then
            tg0=tg0+dtg
            eps1=abs(y2(n)-yk)
            f=y2
        end if
end do
print*,f(1000)
print*,"alfa=",tg0
print*,"dalfa=",dtg
end program solverdiffeq