最近涉及到使用Fortran对微分方程求解，我们知道MATLAB已有内置的函数，比如ode家族，ode15s，对应着不同的求解办法。通过查看odepack的官方文档，我尝试使用了dlsode求解刚性和非刚性常微分方程组。
首先是github网址：https://github.com/jacobwilliams/odepack
具体使用办法：
1.我使用的是vs2022，比较简单的用法就是把，src文件夹所有的文件复制到和项目一个文件夹即可，将M_odepack.f90文件放入到项目中，这样就可以用了。
2.在使用前要use m_odepack
3.这里以官方文档中的例子为例：
program dlsode_ex
use m_odepack
implicit none
external fex
external jex
integer,parameter :: dp=kind(0.0d0)
real(kind=dp),dimension(3) :: atol,y
integer :: iopt,iout,istate,itask,itol,liw,lrw,mf,neq
integer,dimension(23) :: iwork
real(kind=dp) :: rtol,t,tout
real(kind=dp),dimension(58) :: rwork
neq = 3
y(1) = 1.D0
y(2) = 0.D0
y(3) = 0.D0
t = 0.D0
tout = .4D0
itol = 2
rtol = 1.D-4
atol(1) = 1.D-6
atol(2) = 1.D-10
atol(3) = 1.D-6
itask = 1
istate = 1
iopt = 0
lrw = 58
liw = 23
mf = 21
do iout = 1,12
call dlsode(fex,[neq],y,t,tout,itol,[rtol],atol,itask,istate,iopt, &
& rwork,lrw,iwork,liw,jex,mf)
write (6,99010) t,y(1),y(2),y(3)
99010 format (' At t =',d12.4,' y =',3D14.6)
if ( istate<0 ) then
write (6,99020) istate
99020 format (///' Error halt.. ISTATE =',i3)
stop 1
else
tout = tout*10.D0
endif
enddo
write (6,99030) iwork(11),iwork(12),iwork(13)
99030 format (/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4)
end program dlsode_ex
subroutine fex(Neq,T,Y,Ydot)
implicit none
integer,parameter :: dp=kind(0.0d0)
integer :: Neq
real(kind=dp) :: T
real(kind=dp),intent(in),dimension(3) :: Y
real(kind=dp),intent(inout),dimension(3) :: Ydot
Ydot(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)
Ydot(3) = 3.D7*Y(2)*Y(2)
Ydot(2) = -Ydot(1) - Ydot(3)
end subroutine fex
subroutine jex(Neq,T,Y,Ml,Mu,Pd,Nrpd)
implicit none
integer,parameter :: dp=kind(0.0d0)
integer :: Neq
real(kind=dp) :: T
real(kind=dp),intent(in),dimension(3) :: Y
integer :: Ml
integer :: Mu
real(kind=dp),intent(inout),dimension(Nrpd,3) :: Pd
integer,intent(in) :: Nrpd
Pd(1,1) = -.04D0
Pd(1,2) = 1.D4*Y(3)
Pd(1,3) = 1.D4*Y(2)
Pd(2,1) = .04D0
Pd(2,3) = -Pd(1,3)
Pd(3,2) = 6.D7*Y(2)
Pd(2,2) = -Pd(1,2) - Pd(3,2)
end subroutine jex一些变量意义具体看文档说明：https://jacobwilliams.github.io/odepack/proc/dlsode.html
其中，假设n是方程个数，
y：是初值，数组，y(n)
atol:每个方程的绝对误差，数组，atol(n)
t:输入的初始点，tout是下一个点。
mf:是求解方法，其中如果等于21,24需要使用者自己提供雅各比矩阵，如示例代码中jex函数中那样，如果等于10,22,25则不需要自己写，但是jex函数还是需要定义，就是函数框架，函数名，变量声明就可。
fex函数：写的就是你的微分方程组
另外，
rwork，iwork也是两个一维数组，大小如图所示。
以及，
lrw = 22 + 9*NEQ + NEQ**2
liw = 20 + NEQ
整体使用的逻辑就是先设置t值，然后设置循环，tout不断累加，下次循环就使用上次计算得到的新y值以及tout进行迭代计算。
istate是用于输入和输出以指定计算状态的索引，要注意的是如果istate选择2，或3需要在第一次循环中等于1，初始化，到了第二次循环开始才赋值为2或3。