Why not just extend the `lsoda_class`

derived type? Or you can also just specify a function which is contained in another function. Here are the two ways I would do it:

```
module example
use odepack_mod, only: lsoda_class, dp
implicit none
! some data to be passed
type :: MyData
integer :: n
real(dp), allocatable :: arr(:)
end type
! customize the lsoda class
type, extends(lsoda_class) :: lsoda_custom
type(MyData), pointer :: dat => null()
end type
contains
subroutine test_1()
type(lsoda_custom) :: ls
integer :: neq, itask, istate
real(dp) :: y(2), t, tout, rtol, atol(1)
type(MyData), target :: dat
! custom data
dat%n = 3
allocate(dat%arr(dat%n))
dat%arr(:) = 1.0_dp
neq = 2
call ls%initialize(rhs_1, neq, istate=istate)
if (istate < 0) then
error stop '"test" failed'
endif
ls%dat => dat
y(:) = [5.0_dp, 0.8_dp]
t = 0.0_dp
tout = 10.0_dp
rtol = 1.0e-8_dp
atol = 1.0e-8_dp
itask = 1
istate = 1
call ls%integrate(y, t, tout, rtol, atol, itask, istate)
if (istate < 0) then
error stop '"test" failed'
endif
print*,y
end subroutine
subroutine rhs_1(self, neq, t, y, ydot, ierr)
class(lsoda_class), intent(inout) :: self
integer, intent(in) :: neq
real(dp), intent(in) :: t
real(dp), intent(in) :: y(neq)
real(dp), intent(out) :: ydot(neq)
integer, intent(out) :: ierr
type(MyData), pointer :: dat
select type (self)
class is (lsoda_custom)
dat => self%dat
end select
ydot(1) = y(1)-y(1)*y(2)*sum(dat%arr)
ydot(2) = y(1)*y(2)-y(2)
ierr = 0
end subroutine
subroutine test_2()
type(lsoda_class) :: ls
integer :: neq, itask, istate
real(dp) :: y(2), t, tout, rtol, atol(1)
type(MyData) :: dat
! custom data
dat%n = 3
allocate(dat%arr(dat%n))
dat%arr(:) = 1.0_dp
neq = 2
call ls%initialize(rhs_tmp, neq, istate=istate)
if (istate < 0) then
error stop '"test" failed'
endif
y(:) = [5.0_dp, 0.8_dp]
t = 0.0_dp
tout = 10.0_dp
rtol = 1.0e-8_dp
atol = 1.0e-8_dp
itask = 1
istate = 1
call ls%integrate(y, t, tout, rtol, atol, itask, istate)
if (istate < 0) then
error stop '"test" failed'
endif
print*,y
contains
subroutine rhs_tmp(self_, neq_, t_, y_, ydot_, ierr_)
class(lsoda_class), intent(inout) :: self_
integer, intent(in) :: neq_
real(dp), intent(in) :: t_
real(dp), intent(in) :: y_(neq_)
real(dp), intent(out) :: ydot_(neq_)
integer, intent(out) :: ierr_
call rhs_2(neq_, t_, y_, ydot_, ierr_, dat)
end subroutine
end subroutine
subroutine rhs_2(neq, t, y, ydot, ierr, dat)
integer, intent(in) :: neq
real(dp), intent(in) :: t
real(dp), intent(in) :: y(neq)
real(dp), intent(out) :: ydot(neq)
integer, intent(out) :: ierr
type(MyData), intent(inout) :: dat
ydot(1) = y(1)-y(1)*y(2)*sum(dat%arr)
ydot(2) = y(1)*y(2)-y(2)
ierr = 0
end subroutine
end module
program main
use example, only: test_1, test_2
implicit none
call test_1()
call test_2()
end program
```