Maybe you’d like to give ForSolver a try.

Here’s an example using the Newton method:

##
newton

```
program newton
use kinds
use forsolver
implicit none
type(nlsolver) :: nls
real(rk) :: x
call nls%set_options(nl_method='newton', maxit=100, TolFun=1e-6_rk)
call nls%solve(F=F, dFdx=dFdx, x0=0.0_rk, x_sol=x)
contains
function r(x,c1,c2,c3)
real(rk), intent(in) :: x
real(rk), intent(in) :: c1, c2, c3
real(rk) :: r
r = c1 / (c2 + c3*x)
end function
function drdx(x,c1,c2,c3)
real(rk), intent(in) :: x
real(rk), intent(in) :: c1, c2, c3
real(rk) :: drdx
drdx = -c1 * c3 / (c2 + c3 * x)**2
end function
function F(x)
real(rk), intent(in) :: x
real(rk) :: F
real(rk) :: c1, c2, c3
c1 = 10.0_rk
c2 = 12.0_rk
c3 = 30.0_rk
F = r(x,c1,c2,c3) * x**4 + r(x,c1,c2,c3)**2 + r(x,c1,c2,c3)
end function
function dFdx(x)
real(rk), intent(in) :: x
real(rk) :: c1, c2, c3
real(rk) :: dFdx
real(rk) :: dr_dx
c1 = 10.0_rk
c2 = 12.0_rk
c3 = 30.0_rk
dFdx = drdx(x,c1,c2,c3) * x**4 + r(x, c1, c2, c3) * 4.0_rk * x**3&
+ 2.0_rk*drdx(x,c1,c2,c3)*r(x, c1, c2, c3)&
+ drdx(x,c1,c2,c3)
end function
end program
```

##
Results

```
-----------------------------------------------
maxit x0 tol
100 0.00000000 0.1000E-05
-----------------------------------------------
start newton
-----------------------------------------------
it xn F(xn) dF(xn)/dxn
0 0.0000 0.1528E+01 -0.5556E+01
1 0.2750 0.7405E+00 -0.1417E+01
2 0.7975 0.4684E+00 0.1089E+00
3 -3.5056 -0.1631E+02 0.1325E+02
4 -2.2749 -0.4908E+01 0.5772E+01
5 -1.4246 -0.1559E+01 0.2344E+01
6 -0.7592 -0.3751E+00 0.2978E+01
7 -0.6332 0.3838E+00 0.1186E+02
8 -0.6656 0.7384E-01 0.7689E+01
9 -0.6752 0.4195E-02 0.6838E+01
10 -0.6758 0.1534E-04 0.6788E+01
11 -0.6758 0.2069E-09 0.6788E+01
-----------------------------------------------
end newton
-----------------------------------------------
x_sol = -0.67580449156891309
```

And here’s the same example using the complex step Newton method:

##
newton-quasi-cs

```
program newton_quasi_cs
use kinds
use forsolver
implicit none
type(nlsolver) :: nls
complex(rk) :: x
call nls%set_options(nl_method='newton-quasi-cs', maxit=100, TolFun=1e-6_rk, cs_tol=tiny(0.0_rk))
call nls%solve(F=F, x0=(0.0_rk, 0.0_rk), x_sol=x)
contains
function r(x,c1,c2,c3)
complex(rk), intent(in) :: x
complex(rk), intent(in) :: c1, c2, c3
complex(rk) :: r
r = c1 / (c2 + c3*x)
end function
function F(x)
complex(rk), intent(in) :: x
complex(rk) :: F
complex(rk) :: c1, c2, c3
c1 = (10.0_rk, 0.0_rk)
c2 = (12.0_rk, 0.0_rk)
c3 = (30.0_rk, 0.0_rk)
F = r(x,c1,c2,c3) * x**4 + r(x,c1,c2,c3)**2 + r(x,c1,c2,c3)
end function
end program
```

##
Results

```
-----------------------------------------------
maxit x0 tol
100 0.00000000 0.1000E-05
-----------------------------------------------
start newton
-----------------------------------------------
it xn F(xn) dF(xn)/dxn
0 0.0000 0.1528E+01 -0.5556E+01
1 0.2750 0.7405E+00 -0.1417E+01
2 0.7975 0.4684E+00 0.1089E+00
3 -3.5056 -0.1631E+02 0.1325E+02
4 -2.2749 -0.4908E+01 0.5772E+01
5 -1.4246 -0.1559E+01 0.2344E+01
6 -0.7592 -0.3751E+00 0.2978E+01
7 -0.6332 0.3838E+00 0.1186E+02
8 -0.6656 0.7384E-01 0.7689E+01
9 -0.6752 0.4195E-02 0.6838E+01
10 -0.6758 0.1534E-04 0.6788E+01
11 -0.6758 0.2069E-09 0.6788E+01
-----------------------------------------------
end newton
-----------------------------------------------
x_sol = -0.67580449156891309
```