SUBROUTINE matgenph (u, z, w, n, rstart, rend, nstpl, xksq)
USE FileUnits_Module
USE tlogd_Module
USE Storage_Module
USE engpro_Module
USE Masses_Module
USE Numbers_Module
USE dip_Module
IMPLICIT NONE
LOGICAL nsave, wrpotmat
INTEGER n, nstpl, nstep, istep
REAL(Kind=WP_Kind) rstart, rend, r, hstep, d1, d2, d4, drnow, wavel, h, d3
REAL(Kind=WP_Kind) u(1), z(n,n), w(n,n), xksq(n)
EXTERNAL potmatph, smxinv
!            I N T R I N S I C    F U N C T I O N S
intrinsic sqrt
DATA hstep/1.d-4/

nsave=numder
numder=.true.
wrpotmat=.false.

!-----------------------------------------------------------------------
! Read in energy independent information if possible.
!-----------------------------------------------------------------------
wavel=2.d0*pi/sqrt(usys2*engmax)
IF(wavel>2.d0) wavel=2.d0
nstep=nstpl*(rend-rstart)/wavel+.9d0
IF(nstep==0) nstep=1
drnow=(rend-rstart)/nstep
h=drnow*0.5d0
d1=h*h/3.d0
d2=2.d0*d1
d4=-d1/16.d0
d3=2.d0*d2
IF(iwrindep)THEN
   WRITE(frst_unit)wavel,nstep,drnow,h,d1,d2,d4
   IF(lphoto)THEN
      WRITE(fph_unit)wavel,nstep,drnow,h,d1,d2,d3
   ENDIF
ENDIF

!-----------------------------------------------------------------------
! Calculate the interaction matrix.
!-----------------------------------------------------------------------
r=rstart
IF(wrpotmat)THEN
    CALL potmatph (n, r, u, hstep, -1)
ELSE
    WRITE(frst_unit)n,r,hstep
ENDIF
!-----------------------------------------------------------------------
! start integrating the differential equations.
!-----------------------------------------------------------------------
DO istep = 1, nstep
   r=rstart+h*(2*istep-1)
   IF(wrpotmat)THEN
      CALL potmatph (n,r,u,hstep,-1)
   ELSE
      WRITE(frst_unit)n,r,hstep
   ENDIF
   r=rstart+h*(2*istep)
   IF(wrpotmat)THEN
      CALL potmatph (n, r, u, hstep, -1)
   ELSE
      WRITE(frst_unit)n,r,hstep
   ENDIF
   IF(istep==nstep) d2=d1
ENDDO

numder=nsave

RETURN
END