!********************************************************************!

! The definition of the function: f(z) = exp(3*z) + 2*z*cos(z) - 1
subroutine func(z, f, p)

implicit none

integer, parameter :: prec = 8   ! prec = 4 for single precision and prec = 8 for double precision
complex(prec) :: z, f            ! z - value of independent complex variable, f - the function value at z
integer :: p                     ! p - additional (dummy) parameter

f = exp(3.0*z) + 2.0 * z * cos(z) - 1.0

end subroutine

!********************************************************************!

! The definition of the function derivative: df(z)/dz = 3*exp(3*z) + 2*cos(z) - 2*z*sin(z)
subroutine dfunc(z, df, p)

implicit none

integer, parameter :: prec = 8   ! prec = 4 for single precision and prec = 8 for double precision
complex(prec) :: z, df           ! z - value of independent complex variable, df - the derivative value at z
integer :: p                     ! p - additional (dummy) parameter

df = 3.0 * exp(3.0*z) + 2.0 * cos(z) - 2.0 * z * sin(z)

end subroutine

!********************************************************************!

program test

implicit none

integer, parameter :: pp = 8    ! pp = 4 for 32 bit machines and pp = 8 for 64 bit machines
integer, parameter :: prec = 8  ! prec = 4 for single precision and prec = 8 for double precision

integer(pp) :: solver           ! holds pointer to the ZerSol object

integer, parameter :: max_n_zeros = 32         ! specify maximum number of wanted zeros
complex(prec), dimension(max_n_zeros) :: Z, V  ! arrays for zeros and values of the function

integer :: n_zeros                     ! number of zeros
integer :: status                      ! the search status
integer :: n                           ! a counter
integer :: data = 0                    ! additional (dummy) parameter for passing to f(z) and df(z)/dz
real(prec) :: xmin = -2.0, xmax = 2.0  ! x boundaries of the rectangular region
real(prec) :: ymin = -2.0, ymax = 3.0  ! y boundaries of the rectangular region
character(1) :: file = ''              ! name of a file to print the solver info: dimension = len(trim(file)) + 1

external func, dfunc                   ! external subroutines for f(z) and df(z)/dz evaluation

integer(pp) :: f, df                   ! addresses of f(z) and df(z)/dz subroutines

f = loc(func)    ! address of f(z) subroutine
df = loc(dfunc)  ! address of df(z)/dz subroutine, set df = 0 if df(z)/dz is unavailable (a numerical approximation will be used instead)

! pass the function and the region parameters:
call create_solver(solver, f, df, data, xmin, xmax, ymin, ymax)

! the solver sets Z, V, n_zeros and returns 0 if successful:
call find_zeros(solver, max_n_zeros, Z, V, n_zeros, status)

! the solver prints information about the search status on the screen since file is an empty string
call print_status(solver, file, len(trim(file)))

! free the memory allocated by the solver:
call free_solver(solver)

end program

!********************************************************************!