subroutine fcn(t, y, f_bdf)
use work
! Calculeaza functia din BDF
implicit real(8) (a-h,o-z)
real(8) y(n), f_bdf(n), ff(n)

! y =xa(1,:)
! Membrii 2:
!      Fcn(1) = 1.0d0/y(1) - x**2 - 2.0d0/x**3
!      Fcn(2) = y(1)/y(2)**2 - 1.0d0/x - 0.5d0/SQRT(x**3)

ff(1) =1.0d0/y(1) - t**2 - 2.0d0/t**3
ff(2) =y(1)/y(2)**2 - 1.0d0/t - 0.5d0/dSQRT(t**3)

!f_bdf(:) = c(5,0)*y(:) +c(5,1)*xa(0,:) +c(5,2)*xa(-1,:) +c(5,3)*xa(-2,:) +c(5,4)*xa(-3,:) +c(5,5)*xa(-4,:) -h*ff(:)
f_bdf =-h*ff(:) +c(k,0)*y(:)

do i =1, k
f_bdf =f_bdf +c(k, i)*xa(-i+1,:)
enddo


END subroutine fcn

