VMAT_lab3_24

program lab3
  interface
    subroutine dif(T, Y, YP)
    real T, Y(2), YP(2)
    end
  end interface
  integer, parameter :: NEQN = 2
  integer ::  IFLAG, IWORK(6)
  real :: Y(NEQN),YP(NEQN), T, TOUT, RELERR, ABSERR, WORK(3+6*NEQN + 1 ), Y2(NEQN), YP2(NEQN), h, step
  Y(1) = 5.0
  Y(2) = -1.0
  T = 0.0
  RELERR = 0.1E-03
  ABSERR = 0.0
  IFLAG = 1
  HPRINT = 0.05
  TOUT = T
  open(1, file="rkf45.txt")
  do while(TOUT<1.02)
    CALL RKF45(dif,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK)
    write(1,*) "Y = ", Y
    write(1,*) "T = ", T
    write(1,*) "--------------------------"
    TOUT = TOUT + HPRINT
  end do
  close(1)
  Y(1) = 5.0
  Y(2) = -1.0
  T = 0.0
  h = 0.0025
  step = 0.0
  open(2, file="rung.txt")
  do while(T<1.0)
    call  dif(T, Y, YP)
    Y2 = Y + (h/6.0) * YP
    T = T + h/6.0
    call dif(T, Y2, YP2)
    Y = Y + h * (3 * YP2 - 2 * YP)
    if (step >= hprint) then
      step = 0
      write(2,*) "Y = ", Y
      write(2,*) "T = ", T
      write(2,*) "--------------------------"
    end if
    T = T + h - h/6.0
    step = step + h
  end do
  close(2)
  write(*,*) "Values imported to rkf45.txt, rung.txt"
end program lab3

subroutine dif(T, Y, YP)
real T, Y(2), YP(2)
YP(1) = -45 * Y(1) + 60 * Y(2) + sin(1 + T)
YP(2) = 70 * Y(1) - 110 * Y(2) + cos(1 - T) + T + 1
end