[数值问题] 大神，参数赋值的问题感觉没有赋到表达时当中

Module M_rk4_ods
  Implicit None
  Integer :: M = 36000
Contains

  Subroutine Solve(Func, T0, Tt, Y0, N)
    Implicit None
    External Func
    Integer :: N
    Real *8 :: Y0(N), Y(N)
    Real *8 :: K1(N), K2(N), K3(N), K4(N)
    Real *8 :: H, Tt, T0, T
    Integer :: I
    H = (Tt-T0)/M
    T = T0
    Y = Y0
    Do I = 1, M
      Call Func(K1, T, Y, N)
      Call Func(K2, T+H/2, Y+H/2*K1, N)
      Call Func(K3, T+H/2, Y+H/2*K2, N)
      Call Func(K4, T+H, Y+H*K3, N)
      Y = Y + (K1+2*K2+2*K3+K4)*H/6
      T = T0 + I*H
      Write (11, 100) T, Y(1:6)
    End Do
    100 Format (1X, F12.6, 6F15.5)
  End Subroutine Solve

  Subroutine Fun1(F, T, Y, N)
    Implicit None
    Integer :: N
    Real *8 :: F(N), Y(N)
    Real *8 :: T, T0, Tt
    Real *8 :: W, R, Fai, Fai0
    Real *8 :: Ex, Ey, Ez, Bx, By, Bz
    Real Q, W0, K
    Complex I
    I = (0, 1)
    W0 = 200
    Q = 100
    Fai0 = 0.811
    W = W0*Sqrt(1+((2*Y(6))/(W0**2))**2)
    R = Y(6)*(1+((W0**2)/(2*Y(6)))**2)
    Fai = Atan((W0**2)/(2*Y(6)))
    Ex = (W0/W)*Exp((-(Y(4)**2+Y(5)**2)/(W**2)))*Exp(I*(T-Y(6)-Fai-Fai0-((Y(4)**2+Y(5)**2)/(2*R))))
    Ey = 0
    Ez = ((2*Y(4)**2)/(W**2*R))*Ex
    Bx = ((4*Y(4)**2*Y(5)**2)/(W**4*R**2))*Ex
    By = ((1-((W0**2)/(2*Y(6)*R))+((8*Y(6)**3*(Y(4)**2+Y(5)**2)*R)/((R-2*Y(6))**2))-((4*Y(4)**4)/(W**4*R**2)))+I*(((4*Y(6))/(W0**2*W**2))+((4*Y(4))/(W**2*R))-((8*Y(6)*(Y(4)**2+Y(5)**2))/(W0**2*W**4))))*Ex
    Bz = -((2*Y(5)**2)/(W**2*R))*Ex
    F(1) = Q*(Ex+(Y(2)/Sqrt(1+Y(2)**2))*Bz-(Y(3)/Sqrt(1+Y(3)**2))*By)
    F(2) = Q*(Ey+(Y(3)/Sqrt(1+Y(3)**2))*Bx-(Y(1)/Sqrt(1+Y(1)**2))*Bz)
    F(3) = Q*(Ez+(Y(1)/Sqrt(1+Y(1)**2))*By-(Y(2)/Sqrt(1+Y(2)**2))*Bx)
    F(4) = Y(1)/Sqrt(1+Y(1)**2)
    F(5) = Y(2)/Sqrt(1+Y(2)**2)
    F(6) = Y(3)/Sqrt(1+Y(3)**2)
  End Subroutine Fun1

End Module M_rk4_ods

Program Main
  Use M_rk4_ods
  Implicit None
  Integer :: N
  Real *8 :: Y(6), Y0(6)
  Real *8 :: T0, Tt
  Real Px0, Pz0, X0, Z0
  Px0 = 6
  Pz0 = 60
  X0 = -1800
  Z0 = -18000
  Open (Unit=11, File='result0.txt')
  N = 6
  T0 = 0.0
  Tt = 7200.0
  Y0 = (/ Px0, 0, Pz0, X0, 0, Z0 /)
  Call Solve(Fun1, T0, Tt, Y0, N)
End Program Main