Pari/GP, 505 299 196 188 bytes

Saved 317 bytes thanks to @ceilingcat

Golfed version. Try it online!

f(A,m)={local(n,e,s,o);e=matrix(n=#A);for(i=1,n,for(j=1,n,s=matrix(o=n-1);for(r=1,o,for(c=1,o,s[r,c]=A[r+(r>=i),c+(c>=j)]));if((e[j,i]=(-1)^(i+j)*matdet(s)%m)<0,e[j,i]+=m)));e/matdet(A)%m}

Ungolfed version. Try it online!

/* Returns the inverse of matrix A modulo m, or error if not invertible
   Uses the adjugate matrix method:
   A^(-1) ≡ adj(A) * det(A)^(-1) (mod m) */
matinversemod(A, m) = {
    local(n, det, detinv, adj);
    
    if (type(A) != "t_MAT",
        error("Input must be a matrix")
    );
    
    n = matsize(A)[1];
    if (n != matsize(A)[2],
        error("Matrix must be square")
    );
    
    /* Calculate determinant mod m */
    det = matdet(A) % m;
    
    /* Check if determinant is invertible mod m */
    if (gcd(det, m) != 1,
        error("Matrix is not invertible modulo " m)
    );
    
    /* Calculate inverse of determinant mod m */
    detinv = lift(Mod(1/det, m));
    
    /* Calculate adjugate matrix */
    adj = matrix(n,n);
    for(i=1, n,
        for(j=1, n,
            /* (-1)^(i+j) * determinant of minor */
            adj[j,i] = (-1)^(i+j) * matdet(matminor(A,i,j)) % m;
            if (adj[j,i] < 0, adj[j,i] += m)
        )
    );
    
    /* Multiply adjugate by inverse determinant mod m */
    return(lift(Mod(adj * detinv, m)));
}

/* Helper function to get minor matrix by removing row i and column j */
matminor(A,i,j) = {
    local(n, minor);
    n = matsize(A)[1];
    minor = matrix(n-1,n-1);
    
    /* Copy elements, skipping row i and column j */
    for(r=1, n-1,
        for(c=1, n-1,
            minor[r,c] = A[if(r<i,r,r+1), if(c<j,c,c+1)];
        )
    );
    return(minor);
}

Charcoal, 41 34 bytes

ＩＥη⊟ΦＥＸθη÷λＸθ…⁰η⬤λ⁼⁼ιξ﹪ΣＥλ×π§§ζρξθ

Try it online! Link is to verbose version of code. Takes input as \$ m, n, M \$ where \$ n \$ is the size of \$ M \$, and does not reduce its output modulo \$ m \$ (can be done at a cost of 2 bytes). Slow, so times out for the large test cases. (Will crash if given a singular matrix.) Explanation:

  η                                 Input `n`
 Ｅ                                  Map over implicit range
       θ                            Input `m`
      Ｘ                             Raised to power
        η                           Input `n`
     Ｅ                              Map over implicit range
          λ                         Current value
         ÷                          Vectorised integer divide
            θ                       Input `m`
           Ｘ                        Vectorised raise to power
             …                      Range from
              ⁰                     Literal integer `0`
               η                    To input `n`
    Φ                               Filtered where
                 λ                  Current row
                ⬤                   All indices satisfy
                    ι               Outermost value
                   ⁼                Equals
                     ξ              Inner index
                  ⁼                 Equals
                         λ          Current row
                        Ｅ           Map over values
                           π        Current value
                          ×         Multiplied by
                              ζ     Input `M`
                             §      Indexed by
                               ρ    Innermost index
                            §       Indexed by
                                ξ   Inner index
                       Σ            Take the sum
                      ﹪             Modulo
                                 θ  Input `m`
   ⊟                                Take the assumed match
Ｉ                                   Cast to string
                                    Implicitly print

Python 3 + SymPy, 33 bytes

from sympy import*
Matrix.inv_mod

Try it online!

SymPy's Matrix class has a method for modular inverse.

J, ¹⁸ 16 bytes

(]%1+.]^5 p:[)%.

Try it online!

Resolves p/q mod n element-wise (instead of using det(M) to resolve the modular inverse globally).

Abuses GCD of rational numbers to extract 1/q from p/q.

How it works

(]%1+.]^5 p:[)%.    NB. left arg = modulo, right arg = matrix
(            )%.    NB. bind inv(matrix) as new right arg
        5 p:[       NB. phi(modulo)
      ]^            NB. inv(matrix)^phi(modulo) element-wise
   1+.              NB. GCD with 1; GCD(1, p/q) = 1/q
 ]%                 NB. Divide inv(matrix) by the above element-wise

J, 18 bytes

%.@]*-/ .*@]^5 p:[

Try it online!

A dyadic tacit function that takes modulo (left arg) and the matrix (right arg), and gives possibly very large-valued modular inverse of the matrix. To reduce the range, prepend [| at the start of the function.

How it works: the math

A simple mathematical way to calculate the modular inverse of a matrix is the following:

$$ \begin{align} M^{-1} \text{ mod }n &= \text{cofactor}(M) \times \bigl((\det M)^{-1} \text{ mod }n \bigr) \\ &= M^{-1} \times \det M \times \bigl((\det M)^{-1} \text{ mod }n \bigr) \end{align} $$

If the matrix \$M\$ is invertible modulo \$n\$, we know that \$(\det M)^{-1} \text{ mod }n\$ exists, and it can be found using Euler's theorem:

$$ (\det M)^{-1} \equiv (\det M)^{\varphi(n)-1} \text{ mod }n $$

Then we can simplify the original equation to

$$ \begin{align} M^{-1} \text{ mod }n &= M^{-1} \times \det M \times \bigl((\det M)^{\varphi(n)-1} \text{ mod }n \bigr) \\ &\equiv M^{-1} \times (\det M)^{\varphi(n)} \mod{n} \end{align} $$

And now the fun fact: J has built-ins for matrix inverse, matrix determinant, and Euler's totient function. And it uses built-in rational numbers when calculating the matrix inverse!

How it works: the code

%.@]*-/ .*@]^5 p:[    NB. left arg = modulo, right arg = matrix
             5 p:[    NB. totient(modulo)
     -/ .*@]          NB. det(matrix)
            ^         NB. det(matrix) ^ totient(modulo)
%.@]                  NB. inv(matrix)
    *                 NB. inv(matrix) * det(matrix) ^ totient(modulo)

MATL, (25?)  31 29  26 bytes

My first MATL answer

-5 bytes & a bug-fix (+2) thanks to Luis Mendo!

The trailing . may be unnecessary - it is if there is only ever a single inverse of M with elements modulo m.

:inZ^!"&G@[]eY*w\tZyXy=?@.

A full program which prints the elements in row major order separated by newlines.

Try it online! - Too slow for any of the given test cases.

Quite possibly not the best approach for MATL.

How?

:inZ^!"&G@[]eY*w\tZyXy=?@. - expects inputs m and M
:                          - range (m) -> [1,2,...,m]
 i                         - push input (M)
  n                        - number of elements
   Z^                      - ([1,2,...,m]) Cartesian power (#elements(M))
     !                     - transpose
      "                    - for each column, C:
       &G                  -   push both inputs
         @                 -   push C
          []               -   push an empty array (to make e work as below)
            e              -   reshape (C) to square matrix of side ceil(#elements(C)^0.5)
             Y*            -   (reshaped C) matrix multiplication (copy of M)
               w           -   swap top two stack entries
                \          -   (multiplication result) modulo (copy of m)
                 t         -   duplicate top of stack
                  Zy       -   size
                    Xy     -   (size by size) identity matrix
                      =    -   equal -> logical matrix
                       ?   -   if all are truthy:
                        @  -     push C
                         . -     break
                           - implicit print of stack (the valid C)

Jelly, (21?) 22 bytes

The trailing Ṫ may be unnecessary - it is if there is only ever a single inverse of M with elements modulo m.

Ḷṗ⁹L²¤ṁ€⁹æ×%³L⁼þ`$ƑɗƇṪ

A full program printing the result.

Try it online! - Too slow for any of the given test cases (the 35 case took ~20 minutes locally).

11 bytes (but floating point output):

Using Bubler's observation (go upvote!) that raising the determinant to Euler's totient is enough to remove the determinant's denominators:

æ*-×ÆḊ*ÆṪ}ɗ

However, unlike in J, the inversion of \$M\$ in Jelly gives floats so we no longer get an integer matrix as output.

Try it online!

R, 68 bytes

function(M,m,n,A=M){while(any(A%*%M%%m!=diag(n)))A[]=rpois(n^2,9)
A}

Try it online!

Strikingly slow. Will most likely time out for all test cases on TIO, but is guaranteed to give an answer eventually.

Works by rejection sampling: generates random matrices A, with each value taken from a \$Poisson(9)\$ distribution, until a solution is found.

Note that to get A of the correct dimensions, it is 6 bytes shorter to initialize it as A=M and then replace all the values with A[]=rpois(n^2,9) than to create it directly with A=matrix(rpois(n^2,9),n).

R, 97 83 bytes

function(M,m,d){while(any(M%*%(x=matrix(T%/%m^(1:d^2-1),d))%%m-diag(d)))T=T+1;x%%m}

Try it online!

Pretty slow. Takes the dimension of the matrix as input. The previous version using a for loop is a bit faster.

Thanks to Robin Ryder for -14 bytes.

Explanation:

We iterate over every number between \$1\$ and \$m^{d^2}\$, converting each to its base-\$m\$ digits (with leading zeros), reshaping those digits into a matrix of the appropriate size, and testing to see if it's the inverse of \$M\$ modulo \$m\$.

I wanted to attempt the whole series in SNOBOL but I'm not sure I will be able to implement matrix multiplication in SNOBOL in time for it to be a valid submission...

Java 8, 270 261 bytes

M->m->{int l=M.length,R[][]=new int[l][l],T[][]=new int[l][l],d=0,s=l,r,c,k;for(;d!=1|s!=0;){for(r=l*l;r-->0;R[r/l][r%l]=d*=Math.random())d=m;for(d=1,s=r=l;r-->0;d*=T[r][r]%m)for(c=l;c-->0;s-=T[r][c]%m)for(T[r][c]=k=0;k<l;)T[r][c]+=M[r][k]*R[k++][c];}return R;}

-9 bytes thanks to @ceilingcat.

Keeps trying random matrices (including duplicates) until it find the correct one, so times out for most test cases. I tried adding a cache so it tries random matrices without duplicates, but then it still times out for the same test cases.

Try it online (only contains the test cases m=35; M=[[24,14],[48,45]] and m=5; M=[[15,13],[21,13]]).

Explanation:

M->m->{                    // Method with int-matrix & int parameters and int-matrix return
  int l=M.length,          //  Dimension of the input-matrix
      R[][]=new int[l][l], //  Result-matrix of that same size
      T[][]=new int[l][l], //  Temp-matrix of that same size
      d=0,                 //  Flag for the diagonal
      s=l,                 //  Flag for the decreasing sum
      r,c,k;               //  Index integers
  for(;d!=1                //  Continue looping as long as the diagonal flag isn't 1 yet
       |s!=0;){            //  nor the decreasing sum flag isn't 0 yet:
    for(r=l*l;r-->0;       //   Loop over all cells:
      R[r/l][r%l]=         //     Set the current cell in matrix `R`:
        d*=Math.random())d=m;
                           //      To a random value in the range [0,m)
    for(d=1,               //   Reset the diagonal flag to 1
        s=r=l;             //   Reset the decreasing sum flag to `l`
        r-->0              //   Loop over the rows:
        ;                  //     After every iteration:
         d*=               //      Multiply the diagonal flag by:
            T[r][r]        //       The value in the `r,r`'th cell of matrix `T`
                   %m)     //       Modulo the input `m`
      for(c=l;c-->0        //    Inner loop over the columns:
          ;                //      After every iteration:
           s-=             //       Decrease the decreasing sum flag by:
              T[r][c]      //        The value in the `r,c`'th cell of matrix `T`
                     %m)   //        Modulo the input `m`
        for(T[r][c]=k=0;   //     Reset the `r,c`'th cell of matrix `T` to 0
            k<l;)          //     Inner loop `k` in the range [0, length):
          T[r][c]+=        //      Increase the `r,c`'th cell of matrix `T` by:
            M[r][k]        //       The `r,k`'th cell of matrix `M`
            *R[k++][c];}   //       Multiplied by the `k,c`'th cell of matrix `R`
  return R;}               //  And if the loops are done: return matrix `R` as result

Magma, 34 bytes

func<m,M|Matrix(Integers(m),M)^-1>

No TIO for magma, though you can try it on http://magma.maths.usyd.edu.au/calc/

WIN+APL, 114 bytes

Prompts for matrix followed by modulus.

m←r←⎕⋄z←r[1;1]⋄⍎∊(¯1+1↑⍴r)⍴⊂'z←z×1 1↑r←(1 1↓r)-((1↓r[;1])∘.×1↓r[1;])÷r[1;1]⋄'⋄⌊.5+n|((1=n|z×⍳n)/⍳n←⎕)×(z←⌊.5+z)×⌹m

Try it online! Courtesy of Dyalog Classic

SageMath, ⁴⁸ 33 bytes

Saved 15 bytes thanks to ovs!!!

lambda m,M:~Matrix(Integers(m),M)

Nothing on TIO for SageMath unfortunately.

Modular inverse of a matrix M (input as a Python list of lists) mod m.

Sledgehammer, 6 bytes

⠑⡿⡆⠱⣁⣭

Decompresses into this Wolfram Language function:

Inverse[#2, Modulus -> #1]

Try it online!

Jelly, 25 bytes

ÆḊ×Ɱ⁹%ỊTḢ×ZÆḊ-Ƥ$-ƤNÐe⁺€Zʋ

Try it online!

A dyadic link taking the matrix as its left argument and the modulus as its right. Returns a matrix. Append a % to get it within the range 0, m

R, 128 bytes

function(x,m,n)t(round(which((1:m*det(x))%%m<1.5)[1]*outer(1:n,1:n,Vectorize(function(a,b)det(x[-a,-b,drop=F])*(-1)^(a+b))))%%m)

Try it online!

A function taking three arguments, x = the matrix, m = the modulus and n the number of rows of x. Returns a matrix. Uses the same method as my Jelly answer.

JavaScript (ES6),  209  206 bytes

Takes input as (modulo)(matrix).

This transposes the matrix of cofactors (resulting in the adjugate) and multiply it by the inverse of the determinant of \$M\$ modulo \$m\$.

m=>M=>M.map((r,y)=>r.map((_,x)=>((g=k=>(++k*D(M)%m+m)%m-1?g(k):x+y&1?-k:k)``*D(h(M,x).map(r=>h(r,y)))%m+m)%m),h=(a,n)=>a.filter(_=>n--),D=M=>+M||M.reduce((s,[v],i)=>s+(i&1?-v:v)*D(h(M,i).map(r=>h(r,0))),0))

Try it online!

Commented

Helper function \$h\$

The function \$h\$ removes the \$n\$-th entry from the array \$a[\:]\$.

h = (a, n) =>                // a[] = array, n = index
  a.filter(_ => n--)         // keep all but the n-th entry

Helper function \$D\$

The function \$D\$ computes the determinant of the matrix \$M\$.

D = M =>                     // M[] = input matrix
  +M ||                      // if M[] is 1x1, stop recursion and return its unique value
  M.reduce((s, [v], i) =>    // otherwise, for each value v at (0, i):
    s +                      //   add to the sum
    (i & 1 ? - v : v) *      //   either v or -v depending on the parity of i
    D(                       //   multiplied by the result of a recursive call with:
      h(M, i)                //     M[] without the i-th row
      .map(r => h(r, 0))     //     and without the first column
    ),                       //   end of recursive call
    0                        //   start with s = 0
  )                          // end of reduce()

Main function

m => M =>                    // m = modulo, M[] = matrix
  M.map((r, y) =>            // for each position y:
    r.map((_, x) =>          //   for each position x:
      (                      //
        ( g = k =>           //     g is a recursive function taking a counter k
            ( ++k *          //       increment k and multiply it
              D(M)           //       by the determinant of M
              % m + m        //
            ) % m - 1 ?      //       if it's not congruent to 1 modulo m:
              g(k)           //         try again until it is
            :                //       else:
              x + y & 1 ? -k //         return either k or -k
                        : k  //         depending on the parity of x+y
        )`` *                //     initial call to g with a zero'ish value
        D(                   //     multiply by the determinant of:
          h(M, x)            //       M[] without the x-th row
          .map(r => h(r, y)) //       and without the y-th column
        ) % m + m            //     return the result modulo m
      ) % m                  //
    )                        //   end of inner map()
  )                          // end of outer map()

Wolfram Language (Mathematica), 23 bytes

¯\_(ツ)_/¯ the answer was in the documentation of Modulus

Inverse[#2,Modulus->#]&

Try it online!