We are going to reuse all of the code from the beginning of Working with numbers. Here we discuss how to code up the Newton-Raphson method for finding the roots of the cubic van der Waals equation of state problem.

```
.....
! NOTE: everything above here is exactly the same in 16-22a.f90 and 16-22b.f90
write(*,*) "V=? (type the number and press <enter>)"
read (*,*) V ! read initial guess after you look at the plot
write (*,*) '# initial guess: V=',V
write (*,*) '# n','x_n','f(x_n)','fp(x_n)'
do i=1,20
f=c3*(V**3) + c2*(V**2) + c1*V + c0
fp = 3.0*c3*(V**2) + 2.0*c2*V + c1 ! df/dV worked out analytically
write (*,*) i,V,f,fp
V = V-f/fp
end do
write (*,*) '#V = ',V, '(dm^3/mol)'
write (*,*) 'molar density = ',(1/V),'mol/dm^3'
write (*,*) 'mass density = ',(mm/V),'g/dm^3'
```

As usual, let's take it line-by-line:

` read (*,*) V ! read initial guess after you look at the plot`

**reads**from the standard input a value and puts it into the variable V

```
write (*,*) '# initial guess: V=',V
write (*,*) '# n','x_n','f(x_n)','fp(x_n)'
```

```
do i=1,20
f=c3*(V**3) + c2*(V**2) + c1*V + c0
fp = 3.0*c3*(V**2) + 2.0*c2*V + c1 ! df/dV worked out analytically
write (*,*) i,V,f,fp
V = V-f/fp
end do
```

*i*from 1 to 20. When the loop variable is an integer, the default is to increment it by 1. Next, we calculate

*f*and

*fp*, as described in Math Chapter G. (fp is the derivative of f). We print out the values, and then in the last line calculate our new value for V. The new value for "V" is then computed as shown in eqn (G.1) in McQuarrie and Simon.

```
write (*,*) '#V = ',V, '(dm^3/mol)'
write (*,*) 'molar density = ',(1/V),'mol/dm^3'
write (*,*) 'mass density = ',(mm/V),'g/dm^3'
```

The smallest root corresponds to the molar volume of the liquid, and the largest root is the molar volume of the vapor. I used the plot that we computed in Working with numbers to make my initial guess for the values of these "zeros" of the vdW equation. Compiling and running the program above (and typing in "0.0" or "0.3" when appropriate, so that the **read** statement we discussed above can "get" the value), here is my output:

**For the liquid**

```
jschrier@ws7:~/pchem1/ps1$ ./a.out
#conditions = 35.00000 142.6900
#vdW eqn: ( 1.000000 ) v^3 + ( -0.3663687 ) v^2 + (
3.8020000E-02 ) v + ( -1.2101767E-03 )
0.0
# initial guess: V= 0.000000
# nx_nf(x_n)fp(x_n)
1 0.000000 -1.2101767E-03 3.8020000E-02
2 3.1830005E-02 -3.3893756E-04 1.7736409E-02
3 5.0939709E-02 -9.1941329E-05 8.4791277E-03
4 6.1782964E-02 -2.3833243E-05 4.2007118E-03
5 6.7456581E-02 -5.6446297E-06 2.2432059E-03
6 6.9972903E-02 -1.0224758E-06 1.4368519E-03
7 7.0684507E-02 -7.8929588E-08 1.2157112E-03
8 7.0749432E-02 -5.8207661E-10 1.1956841E-03
9 7.0749916E-02 0.000000 1.1955388E-03
10 7.0749916E-02 0.000000 1.1955388E-03
11 7.0749916E-02 0.000000 1.1955388E-03
12 7.0749916E-02 0.000000 1.1955388E-03
13 7.0749916E-02 0.000000 1.1955388E-03
14 7.0749916E-02 0.000000 1.1955388E-03
15 7.0749916E-02 0.000000 1.1955388E-03
16 7.0749916E-02 0.000000 1.1955388E-03
17 7.0749916E-02 0.000000 1.1955388E-03
18 7.0749916E-02 0.000000 1.1955388E-03
19 7.0749916E-02 0.000000 1.1955388E-03
20 7.0749916E-02 0.000000 1.1955388E-03
#V = 7.0749916E-02 (dm^3/mol)
molar density = 14.13429 mol/dm^3
mass density = 564.6367 g/dm^3
FORTRAN STOP
```

**For the vapor**

```
jschrier@ws7:~/pchem1/ps1$ ./a.out
#conditions = 35.00000 142.6900
#vdW eqn: ( 1.000000 ) v^3 + ( -0.3663687 ) v^2 + (
3.8020000E-02 ) v + ( -1.2101767E-03 )
0.3
# initial guess: V= 0.3000000
# nx_nf(x_n)fp(x_n)
1 0.3000000 4.2226380E-03 8.8198751E-02
2 0.2521236 1.1134231E-03 4.3978527E-02
3 0.2268062 2.3375126E-04 2.6153758E-02
4 0.2178686 2.4371315E-05 2.0779714E-02
5 0.2166958 3.9441511E-07 2.0110056E-02
6 0.2166762 -4.6566129E-10 2.0098940E-02
7 0.2166762 4.6566129E-10 2.0098954E-02
8 0.2166762 -4.6566129E-10 2.0098940E-02
9 0.2166762 4.6566129E-10 2.0098954E-02
10 0.2166762 -4.6566129E-10 2.0098940E-02
11 0.2166762 4.6566129E-10 2.0098954E-02
12 0.2166762 -4.6566129E-10 2.0098940E-02
13 0.2166762 4.6566129E-10 2.0098954E-02
14 0.2166762 -4.6566129E-10 2.0098940E-02
15 0.2166762 4.6566129E-10 2.0098954E-02
16 0.2166762 -4.6566129E-10 2.0098940E-02
17 0.2166762 4.6566129E-10 2.0098954E-02
18 0.2166762 -4.6566129E-10 2.0098940E-02
19 0.2166762 4.6566129E-10 2.0098954E-02
20 0.2166762 -4.6566129E-10 2.0098940E-02
#V = 0.2166762 (dm^3/mol)
molar density = 4.615182 mol/dm^3
mass density = 184.3673 g/dm^3
FORTRAN STOP
```