I strongly suggest you look at examples in Section 4.17.4 and 4.2.7 respectively to understand how to perform basin analysis for density difference, and how to refine the positions of the minima/maxima located via basin analysis module by topology analysis module.

If what you need is the exact positions of maxima and minima of deformation density rather than (3,-1) and (3,+1) CPs, the following steps are the best and robust:

There are three steps:
(1) Generate generate grid data of deformation density
Boot up Multiwfn and input path of H2O wavefunction file, then input
5  // Grid data calculation
-2  // Obtain deformation property
1  // Electron density
3  // High quality grid
0  // Return to main menu

(2) Perform basin analysis to locate maxima and minima
17  // Basin analysis modlue
1  // Generate basins
2  // Use grid data in memory
-4  // Export
3  // Export position of located maxima and minima to attractors.txt in current folder
-10  // Return to main menu

(3) Topology analysis on deformation density
iu  // Change user-defined function
-3  // The user-defined function will correspond to interpolation function via B-spline algorithm, as mentioned in Section 2.7 of manual
2  // Topology analysis
-11  // Change the function to be analyzed
100  // User-defined function
1  // Search CPs from given starting points
4  // Using starting points from a .txt file
attractors.txt
0  // Return

Now you can examine current CPs after entering option 0, purple and green spheres correspond to minima and maxima.

Then if you want to compare the CP positions with isosurfaces to confirm if the positions are indeed correct, you can then input

-10  // Return to main menu
13  // Process grid data
-2  // Examine isosurface

After properly adjust isovalue, to e.g. 0.09, as shown below, you can clearly see that green and blue isosurfaces enclose purple and green spheres respectively, so the minima and maxima of deformation density are indeed successfully located.

If you also need (3,-1) and (3,+1) CPs, after locating minima and maxima in topology analysis module, you can input
6  // Search CPs from a batch of points within a sphere
-1  // Start the search using each nucleus as sphere center in turn (I suggest choose this option several times until no new CPs can be found)
-9  // Return

The current CPs shown in option 0:

As you can see, basically all CPs have been found. Satisfaction of poincare-hopf relationship is never important for complicated function like deformation density, you just need to focus on the CPs that you are really interested in.

In fact, it is not impossible to locate minima and maxima without making use of the basin analysis module and then refining their positions, however if you directly perform topology analysis for deformation density via option 6 in topology analysis module, positions of some minima and maxima can hardly be located due to numerical reason and the complexity of this function.

Then with the main option "13" I visualize the CPs and the isosurface (see Fig. 1). My question is whether this sequence of steps is correct for this type of analysis since the poincare-hopf relationship is not satisfied.