^{sandbox (deleted)}

Lets define a matrix of 9s as:
$$ N = begin{bmatrix} 9&9&9\9&9&9\9&9&9 end{bmatrix} $$

Lets define an exploding number as a number at position $(x,y)$ that can be decomposed into equal integers between all its adjacent neighbors (including itself) and the absolute value of each portion is greater than 0.

From the previous matrix, lets explode the number at position $(1,1)$ (0 indexed)
$$ N = begin{bmatrix} 9&9&9\9&color{red}9&9\9&9&9 end{bmatrix} $$
$$ N = begin{bmatrix} 9+color{red}1&9+color{red}1&9+color{red}1\9+color{red}1&color{blue}0+color{red}1&9+color{red}1\9+color{red}1&9+color{red}1&9+color{red}1 end{bmatrix} $$

$$ N = begin{bmatrix} 10&10&10\10&color{red}1&10\10&10&10 end{bmatrix} $$

Sometimes, decomposing result into a rational number greater than 1. This is something we need to avoid when exploding numbers. In this cases the remainder will be assigned to the exploded number.

To demonstrate it, lets continue working with our previous matrix. This time we will explode the number at position $(0,0)$

$$ N = begin{bmatrix} color{red}{10}&10&10\10&1&10\10&10&10 end{bmatrix} $$

Here we have 3 neightbors and the number itself. Here the equation is something like $10/4$ which give us 2 for each and 2 as remainder.

$$ N = begin{bmatrix} color{blue}2+color{red}{2}&color{red}{10+2}&10\color{red}{10+2}&color{red}{1+2}&10\10&10&10 end{bmatrix} $$

$$ N = begin{bmatrix} color{red}{4}&12&10\12&3&10\10&10&10 end{bmatrix} $$

As well, sometimes a number wont be big enough to be decomposed in equal parts (where the abs is greater than 0) between his neighbors (|rational number| < 1). In this cases we need to "borrow" from the exploded number in order to maintain the "greater than 0" condition. Lets continue with our previous example and explode the number at position $(1,1)$.

$$ N = begin{bmatrix} 4&12&10\12&color{red}3&10\10&10&10 end{bmatrix} $$

$$ N = begin{bmatrix} 4+color{red}1&12+color{red}1&10+color{red}1\12+color{red}1&color{blue}0+color{red}{1}-color{green}6&10+color{red}1\10+color{red}1&10+color{red}1&10+color{red}1 end{bmatrix} $$
$$ N = begin{bmatrix} 5&13&11\13&color{red}{-5}&11\11&11&11 end{bmatrix} $$

The challenge is, given a list of $(x,y)$ positions and an finite non-empty array of natural numbers, return the exploded form after each number from the positions list has been exploded.

Test cases

Input: initial matrix: [[3, 3, 3], [3, 3, 3], [3, 3, 3]], numbers: [[0,0],[0,1],[0,2]]

Output: [[1, 0, 1], [5, 6, 5], [3, 3, 3]]

Input: Initial matrix: [[9, 8, 7], [8, 9, 7], [8, 7, 9]], numbers: [[0,0],[1,1],[2,2]]

Output: [[4, 11, 8],[11, 5, 10],[9, 10, 4]]

Input: Initial matrix: [[0, 0], [0, 0]], numbers: [[0,0],[0,0],[0,0]]

Output: [[-9, 3],[3, 3]]

Input: Initial Matrix: [[10, 20, 30],[30, 20, 10],[40, 50, 60]], numbers: [[0,2],[2,0],[1,1],[1,0]]

Output: [[21, 38, 13], [9, 12, 21], [21, 71, 64]]

C(GCC) 220 216 214 bytes

#define L(v)for(int v=2;~v--;)

#define P l/C+r<0|l/C+r>=R|l%C+c<0|l%C+c>=C

f(int R,int C,int*m){for(int*i=m+R*C;~*i;){int*M,l=*i+++C**i++,a=0,b;L(r)L(c)P?:++a;M=m+l;b=*M/a;b+=!b;*M-=b*a;L(r)L(c)*(M+r*C+c)+=P?0:b;}}

Run it here

a slightly less golfed version

#define L(v)for(int v=2;~v--;)

#define P l/C+r<0|l/C+r>=R|l%C+c<0|l%C+c>=C

f(int R, int C, int*m) {

    for(int*i=m+R*C;~*i;) {

        int*M,l=*i+++C**i++,a=0,b;

        L(r)

            L(c)

                P?:++a;

        M=m+l;

        b=*M/a;

        b+=!b;

        *M-=b*a;

        L(r)

            L(c)

                *(M+r*C+c)+=P?0:b;

    }

}

The calling code with an example

int main()

{

  int matrix[] = {3,3,3,3,3,3,3,3,3,0,0,0,1,0,2,-1};

  int rows = 3;

  int columns = 3;

  f(rows,columns,matrix);

  for(int r = 0; r < rows; ++r) {

    for(int c = 0; c < columns; ++c) {

      printf("%03d,",matrix[r*columns + c]);

    }

    printf("n");

  }

}

and the output

001,005,003,

000,006,003,

001,005,003,

Clean, 181 167 bytes

import StdEnv;

foldlm(x,y)={{if(d>2)0b+e-if(d>0)0b*n\e<-:l&v<-[0..],let{b=max m.[y,x]n/n;$a b=2+sign a-(a+1)/size b;n= $x l* $y m;d=(v-x)^2+(u-y)^2}}\l<-:m&u<-[0..]}

Try it online!

In the form of a partially-applied function literal.

Expanded (first version):

f // functinon f on {{Int}} and [(Int,Int)]

    = foldl m (x, y) // fold :: (a -> b -> a) a [b] -> a with first argument  {{Int}} (Int,Int) -> {{Int}} giving  {{Int}} [(Int,Int)] -> {{Int}}

        = {                     // an array of

            {                   // arrays of

                if(d > 2) 0 b   // the amount we give to the neighbors

                + e             // plus the current entry

                - if(d > 0) 0 b // minus the amount taken from the target entry

                * n             // times the number of neighbors, if we're on the target

            \                  // for each

                e <-: l         // element of row l

                & v <- [0..]    // and x-index v

                , let           // local definitions:

                    b           // the amount given to the neighbors

                        = max   // we need at least 1 each, so take the largest of

                            m.[y, x] // the target entry

                            n   // or the number of neighbors

                        / n     // divide it by the number of neighbors

                    n           // the number of neighbors

                        = (     // sum of

                            1   // one

                            + s x // if x is at the left edge = 0 else 1

                            + s ( // if x is at the right edge = 0 else 1

                                size l

                                - x 

                                - 1

                            )

                        ) * (   // times the sum of

                            1   // one

                            + s y // if y is at the top edge = 0 else 1

                            + s ( // if y is at the bottom edge = 0 else 1

                                size m

                                - y

                                - 1

                            )

                        )

                    d           // distance from the target point

                        = (v - x)^2

                        + (u - y)^2

            }

        \                      // for each

            l <-: m             // row l in matrix m

            & u <- [0..]        // and y-index u

        }