As depicted here, what you are looking for is the 1st-order 2-dimensional Voronoi diagram under the Manhattan, or L1-metric. This is a quite non-trivial problem (to solve efficiently), fortunately with many existing algorithms and software.
You actually want the subset of the Voronoi diagram that coincides with the discrete grid defined by your matrix coordinates, but it is straightforward to find which points of the diagram happen to have integral coordinates.
The "fixed points" are called sites, seeds, generators, or source points in different contexts. Since the number of sites is small, it is more appropriate to use a computational geometry algorithm, that is, find a vector representation of the line segments the diagram consists of, independently of the size of your matrix (for 10 sites, it should be extremely fast).
A highly recognized library is CGAL, which includes a large number of computational geometry algorithms in C++, and in particular 2D segment Delaunay graphs (Voronoi diagrams). The user manual talks only about Euclidean distance, but there is also an L1 Voronoi diagram demo, so I hope you'll find your way.
Another perspective is a raster representation, using algorithms that sweep the entire plane (your matrix) once, in a particular order similar to wave propagation. This is recommended if the number of sites is large, forming arbitrary shapes on the plane. It is related to the distance transform, the medial axis and the topological skeleton. The latter two are generalizations of the Voronoi diagram, in that they consist of curves rather than line segments in general, regardless of the metric.
A very good solution under this perspective is distance transforms of sampled functions, also giving source code. This is limited to computing the distance transform and to separable distance metrics for instance squared Euclidean or L1 (your case).
To find the medial axis given the distance transform, a good choice is my own medial feature detector, with binary code available. This work also uses an alternative method for distance transform computation that has no limitation on the distance metric. It is slower than the previous one but still linear-time on the size of domain (your matrix). Computation of the medial axis itself is extremely fast (linear-time on the number of points of the medial axis).
As an example, you can have a look at this image, where white points are the sites (actually forming a closed curve on the plane), and red-yellowish curves are the (Euclidean) medial axis. Note that the latter have sub-pixel accuracy, i.e. appear as 1-pixel-wide anti-aliased curves, although the input sites do not. This illustrates that you have the required information (exact position) to find out whether a candidate point on the medial axis has integral coordinates or not. To get an idea of speed, this result was generated in around 300ms.
If you care about efficiency, my advice is to look for existing algorithms or code. If you insist on your own implementation, you will have to read a lot first, so that you have a good understanding of the problem.