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)

! Ecuatia Van der Pol:
! x'' -lambda(1-x**2)*x'+x =A*cos(omega*t)
! 1st order: x =y(1) & u =y(2)
! x' = u
! u' = lambda*(1-x**2)*u -x + A*cos(omega*t)
!
! x(0) =1;	u(0)=0

! Membrii 2:
ff(1) =y(2)
ff(2) =lambda*(1-y(1)**2)*y(2) -y(1) + A*dcos(omega*t)

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

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

