Charcoal, 81 bytes

ＦθＦι⊞υκＵＭθκ≔ＬθηＦυＦ⁺⁺⪪ＥυληＥθ⁺λ×θη⟦×θ⊕η×⊕θ⊖η⟧«≔Ｅκ§υλι¿⁼¹№ι⁰§≔υ§κ⌕ι⁰⁻÷×⊕×ηηη²Σι»Ｉ⪪υη

Try it online! Link is to verbose version of code. Uses zero as the "blank" marker. Explanation:

ＦθＦι⊞υκ

Flatten the input array.

ＵＭθκ

Replace the original array with a range from 0 to n-1.

≔Ｌθη

Also the length of the array is used a lot so capture it in a temporary to save 3 bytes.

Ｆυ

Loop n² times, which is more than enough to track down all the solvable 0s.

Ｆ⁺⁺

Loop over all of the following ranges:

⪪Ｅυλη

the range from 0 to n²-1, split into subranges of length n;

Ｅθ⁺λ×θη

the subranges obtained from the range from 0 to n²-1, but taking every nth element (so effectively the transpose of the above);

⟦×θ⊕η×⊕θ⊖η⟧«

the range from 0 to n²-1 in steps of n+1, which is the main diagonal, and the range from n-1 to n²-n in steps of n-1, which is the main antidiagonal.

≔Ｅκ§υλι

Get the values in the flattened array corresponding to the elements of the current range.

¿⁼¹№ι⁰

Count whether exactly one of them is zero.

§≔υ§κ⌕ι⁰

If so then overwrite that entry in the flattened array...

⁻÷×⊕×ηηη²Σι

... with ½n(n²+1) minus the sum of the (other) elements.

»Ｉ⪪υη

Split the flattened array back into rows and convert the values to strings for implicit print.

Jelly, 25 bytes

ZṚ,⁸;Jị"$€$§FE
²Œ!ṁ€ÇƇ=ÐṀ

A full program taking n and a list-of-lists-formatted representation of the incomplete square which prints the result in the same format.

Try it online! - too slow for TIO's 60s limit
...so, Try a limited space one which only considers the first 150K permutations - three magic squares two of which match at two and three locations.

How?

Unfortunately, even with the ability to deduce the missing numbers one at a time, I believe brute-forcing will be terser, so that is how this works.

ZṚ,⁸;Jị"$€$§FE - Link 1, Is this a magic-square?: list of lists, M
Z              - transpose (M)
 Ṛ             - reverse (together ZṚ rotate 1/4)
  ,⁸           - pair with chain's left argument (M)
          $    - last two links as a monad:
         €     -   for each (m in (MZṚ, M)):
        $      -     last two links as a monad:
     J         -       range of length = [1..n]
       "       -       zip with:
      ị        -         index into - i.e. get the leading diagonal
    ;          -   concatenate (m with it's diagonal)
           §   - sums
            F  - flatten
             E - all equal?

²Œ!ṁ€ÇƇ=ÐṀ - Main Link: integer, N; list of lists, P
²          - square (n)
 Œ!        - all permutations of (implicit range [1..n²])
   ṁ€      - mould each like (P)
      Ƈ    - filter keep those for which:
     Ç     -   call the last Link as a monad  - i.e. keep magic squares
        ÐṀ - keep those which are maximal under:
       =   -   equals (P) (vectorises) - i.e. keep the one which matches at all givens
           - implicit print, which when given a list containing only one item prints that item

R, 169 180 142 135 bytes

Edits: +11 bytes to rotate the magic square back to it's original orientation, -38 bytes by wrapping "replace-only-missing-element" into a function, -7 bytes by various golf obfuscations

function(m,n){while(F%%4|sum(!m)){m[n:1,]=apply(m,1,f<-function(v){if(sum(!v)<2)v[!v]=(n^3+n)/2-sum(v);v})
m[d]=f(m[d<-!0:n])
F=F+1}
m}

Try it online!

Solves rows & the first diagonal, then rotates the matrix anticlockwise (so cols become rows in the opposite order) and repeats, until there are no empty elements left. Outputs the completed magic-square matrix in one of 4 possible rotated forms. Then rotates matrix back to its original orientation.

Commented readable version:

solve=function(m,n){
    t=(n^3+n)/2                         # t = desired total of each row/col/diag
    f=function(v){                      # f = function to check if a vector
        if(sum(!v)==1)v[!v]=t-sum(v);v  # has only 1 missing element, and if so
    }                                   # fill it with t-sum(elements).
    while(F%%4|sum(!m)){                # While rotations are not multiple-of-4, or
                                        # there are still some empty elements of m:
        m[n:1,]=                        # rotate the matrix anticlockwise, while
            apply(m,1,f)                # using f() to fix any rows; then
        d=1:(n+1)==1                    # define diagonal as every (n+1)th element,
        m[d]=f(m[d])                    # and use f() to fix diagonal.
        F=F+1                           # Count rotations so far,
    }                                   # and repeat.
    m                                   # Finally, output m.
}

or awfully-slow R, 124 123 109 105 bytes

Edit: -14 bytes thanks to Xi'an

function(m,n){x=m;`?`=rowSums;while(any(sum(x[0:n<1])!=c(sum(diag(x)),?x,?t(x))))x[!m]=sample(n^2)[-m];x}

Try it online!

Generates random permutations of the missing elements until it finds one where all row, column and diagonal sums are all equal.

JavaScript (ES7),  143 142  140 bytes

Expects (n)(m), where unknown cells in m are filled with 0's.

n=>g=m=>[0,1,2,3].some(d=>m.some((r,i)=>m.map((R,j)=>t^(t-=(v=d?R:r)[x=[j,i,j,n+~j][d]])||(e--,X=x,V=v),e=1,t=n**3+n>>1)&&!e))?g(m,V[X]=t):m

Try it online!

Commented

n =>                          // outer function taking n
  g = m =>                    // inner function taking the matrix m[]
    [0, 1, 2, 3]              // list of directions
    .some(d =>                // for each direction d:
      m.some((r, i) =>        //   for each row r[] at position i in m[]:
        m.map((R, j) =>       //     for each row R[] at position j in m[]:
          t ^ (               //       test whether t is modified:
            t -=              //         subtract from t:
              (v = d ? R : r) //           use v = r[] if d = 0 or v = R[] otherwise
              [x =            //           use:
                [ j,          //             r[j] if d = 0 (rows)
                  i,          //             R[i] if d = 1 (columns)
                  j,          //             R[j] if d = 2 (diagonal)
                  n + ~j      //             R[n - 1 - j] if d = 3 (anti-diagonal)
                ][d]          //
              ]               //
          ) || (              //       if t was not modified:
            e--,              //         decrement e
            X = x,            //         copy x to X
            V = v             //         copy v to V
          ),                  //
          e = 1,              //       start with e = 1
          t = n**3 + n >> 1   //       start with t = n(n²+1)/2
        )                     //     end of map()
        && !e                 //     e = 0 means that there's exactly one cell set
                              //     to zero in this vector
      )                       //   end of inner some()
    ) ?                       // end of outer some(); if truthy:
      g(m, V[X] = t)          //   update V[X] to t and do a recursive call
    :                         // else:
      m                       //   done: return m[]

APL (Dyalog Unicode), 60 bytes

{(⍵,m+.×1+⍺*2)⌹(∘.(×⊢×=)⍨⍵)⍪2×m←(⍪↑c(⌽c))⍪(⊢⍪⍴⍴⍉)⍺/c←∘.=⍨⍳⍺}

Try it online!

Not likely the shortest approach, but anyway here is one with Matrix Divide ⌹, a.k.a. Solve Linear Equation. This works because all the cells are uniquely determined by the horizontal/vertical/diagonal sums when joined with the givens. ⌹ has no problem with over-determined systems, as long as there is a solution (otherwise, it finds the least-squares fit).

A dyadic inline function (dfn) where the left argument is the side length and the right argument is the flattened matrix.

In the case of [4 9 2][3 0 0][0 0 0], the coefficient matrix and constant vector are given as follows:

Coefficients        Constants
-------------------------------
Part 1: Givens
1 0 0 0 0 0 0 0 0   4
0 1 0 0 0 0 0 0 0   9
0 0 1 0 0 0 0 0 0   2
0 0 0 1 0 0 0 0 0   3
0 0 0 0 0 0 0 0 0   0
0 0 0 0 0 0 0 0 0   0
0 0 0 0 0 0 0 0 0   0
0 0 0 0 0 0 0 0 0   0
0 0 0 0 0 0 0 0 0   0

Part 2: Magic Square sums
2 0 0 0 2 0 0 0 2   30  # diagonals
0 0 2 0 2 0 2 0 0   30
2 2 2 0 0 0 0 0 0   30  # rows
0 0 0 2 2 2 0 0 0   30
0 0 0 0 0 0 2 2 2   30
2 0 0 2 0 0 2 0 0   30  # columns
0 2 0 0 2 0 0 2 0   30
0 0 2 0 0 2 0 0 2   30

which is a set of 17 equations for 9 unknowns.

{(⍵,m+.×1+⍺*2)⌹(∘.(×⊢×=)⍨⍵)⍪2×m←(⍪↑c(⌽c))⍪(⊢⍪⍴⍴⍉)⍺/c←∘.=⍨⍳⍺}

m←(⍪↑c(⌽c))⍪(⊢⍪⍴⍴⍉)⍺/c←∘.=⍨⍳⍺  ⍝ Construct the sums part of the coef matrix
                     c←∘.=⍨⍳⍺  ⍝ ⍺ × ⍺ identity matrix
                   ⍺/  ⍝ ⍺ copies of each horizontally, giving the "rows" part
            (  ⍴⍴⍉)    ⍝ Reshape the transpose of above into the original,
                       ⍝ giving the "columns" part
             ⊢⍪        ⍝ Vertically concatenate two parts
m←(⍪↑c(⌽c))⍪  ⍝ Generate the "diagonals" part and vertically prepend to above

(∘.(×⊢×=)⍨⍵)⍪2×m  ⍝ Construct the entire coef matrix
             2×m  ⍝ Use twos so that we can avoid halving the constant
(          )⍪     ⍝ Vertically concatenate with...
 ∘.(×⊢×=)⍨⍵       ⍝ The square diagonal matrix where nonzero entries of ⍵ give
                  ⍝ a 1 at the corresponding position, 0 otherwise

(⍵,m+.×1+⍺*2)  ⍝ Construct the constant vector
       1+⍺*2   ⍝ Square of ⍺ plus 1
   m+.×        ⍝ Matmul with m, which has ⍺ ones on each row,
               ⍝ giving (# of rows of m) copies of ⍺ times above
 ⍵,            ⍝ Prepend ⍵ to above

⌹  ⍝ Solve the linear system of equations; no postprocessing necessary

Wolfram Language (Mathematica), 100 bytes

#/.Solve[Tr/@Flatten[{#,Thread@#,{(d=Diagonal)@#,d@Reverse@#}},1]==Table[(l^3+l)/2,2(l=Tr[1^#])+2]]&

Try it online!

JavaScript, ⁵⁵⁹ 551 bytes

Fast and methodical.

B=Boolean,f=((e,r)=>(v=r*((r**2+1)/2),e.forEach(e=>e.filter(B).length==r-1?e[e.findIndex(e=>!e)]=v-e.reduce((e,f)=>!(e+=f)||e):0),e[0].reduce((f,l,n)=>!(f[0].push(e[n][n])+f[1].push(e[n][r-1-n]))||f,[[],[]]).forEach((f,l)=>{f.filter(B).length==r-1&&(z=f.findIndex(e=>!e),e[z][l?r-1-z:z]=v-f.reduce((e,f)=>!(e+=f)||e))}),e[0].reduce((f,r,l)=>f.forEach((f,r)=>f.push(e[l][r]))||f,new Array(r).fill().map(()=>[])).forEach((f,l)=>f.filter(B).length==r-1?e[f.findIndex(e=>!e)][l]=v-f.reduce((e,f)=>!(e+=f)||e):0),e.flat(2).filter(B).length==r*r?e:f(e,r)));

Live examples:

B=Boolean,f=((e,r)=>(v=r*((r**2+1)/2),e.forEach(e=>e.filter(B).length==r-1?e[e.findIndex(e=>!e)]=v-e.reduce((e,f)=>!(e+=f)||e):0),e[0].reduce((f,l,n)=>!(f[0].push(e[n][n])+f[1].push(e[n][r-1-n]))||f,[[],[]]).forEach((f,l)=>{f.filter(B).length==r-1&&(z=f.findIndex(e=>!e),e[z][l?r-1-z:z]=v-f.reduce((e,f)=>!(e+=f)||e))}),e[0].reduce((f,r,l)=>f.forEach((f,r)=>f.push(e[l][r]))||f,new Array(r).fill().map(()=>[])).forEach((f,l)=>f.filter(B).length==r-1?e[f.findIndex(e=>!e)][l]=v-f.reduce((e,f)=>!(e+=f)||e):0),e.flat(2).filter(B).length==r*r?e:f(e,r)));

console.log(JSON.stringify(f([
  [4, 9, 2],
  [0, 5, 0],
  [0, 0, 0]
], 3)));

console.log(JSON.stringify(f([
  [1, 23, 0, 4, 21],
  [15, 14, 0, 18, 11],
  [0, 0, 0, 0, 0],
  [20, 8, 0, 12, 6],
  [5, 3, 0, 22, 25]
], 5)));

The "un"-golfed version can be seen at this Github repository.

MATL, 36 bytes

`nZ@[]etGg)GXz-yt!hs&ytXdwPXdhsh&-ha

The input is an $ n times n$ matrix, with $0$ for the unknown numbers.

The code keeps generating random $ n times n$ matrices formed by the numbers $1, dots, n^2$ until one such matrix meets the required conditions. This procedure is guaranteed to finish with probability one.

This is a terrible approach, in that:

• running time is random, and unbounded;
• average running time increases as $(n^2)!$ (that's more than exponentially);
• and so it very likely times out in the online interpreter.

... but hey, it's the shortest answer so far!

(Don't) try it online.

See below a sped-up animated GIF of an example run which happened to take about 2 minutes, here compressed to a few seconds.

Explanation

`         % Do...while
  n       %   Number of elements. This implictly takes the input in the first
          %   iteration, or uses the candidate solution from the previous iteration.
          %   Let this number be denoted as N
  Z@      %   Random permutation of integers 1, 2, ..., N
  []e     %   Reshape as a square matrix. This yields a candidate solution
  t       %   Duplicate
  Gg)     %   Push input, convert to logical, index: this produces a column vector
          %   of the entries of the candidate solution that correspond to nonzero
          %   entries in the input matrix
  GXz     %   Push input, take its nonzero elements. Gives a column vector
  -       %   Element-wise difference (*). This will be all zeros for a valid
          %   solution
  y       %   Duplicate second-top object from the stack, that is, the candidate
          %   solution
  t!      %   Duplicate, transpose
  h       %   Concatenate horizontally
  s       %   Sum of columns. This also gives the sum of rows, thanks to the
          %   concatenated, transposed copy. The result is a two-element row
          %   vector (**)
  &y      %   Duplicate third-top object from the stack: the candidate solution
  tXd     %   Duplicate, extract diagonal as a column vector
  wPXd    %   Swap, flip vertically, extract diagonal. This gives the anti-diagonal
          %   as a column vector
  h       %   Concatenate horizontally
  s       %   Sum of each column. This gives the sum of the diagonal and that  
          %   of the anti-diagonal
  h       %   Concatenate horizontally with (**)
  &-      %   Matrix of all element-wise differences. This will be a matrix of
          %   zeros for a valid solution (***)
  h       %   Concatenate (*) and (***) horizontally. Since sizes do not match,
          %   both (*) and (***) are first linearized to row vectors, and the
          %   result is a row vector
  a       %   Any. This gives true if any element is non-zero
          % End (implicit). A new iteration is run if the top of the stack is true
          % Display (implicit). The candidate solution from the last iteration is
          % the valid solution

05AB1E, 43 41 30 bytes

-2 bytes by replacing Dgt with ¹ to get first input back

-11 bytes thanks to Kevin Cruijssen!

nLœʒ¹ôD©ø®Å®Å/)O˜Ë}ʒøε¬_sË~}P

Try it online! Takes input as (n, flattened square), where zeros represent blanks, like

3
[4,9,2,3,0,0,0,0,0]

Works by generating all permutations of the numbers from 1 to n², filtering to only keep those that are magic squares, and then iterating through and printing all that match the partial input (by input constraints, there will always only be one match). Because of this brute-force approach, it's already very slow for 3x3 magic squares and I doubt 5x5 would terminate. This is my first 05AB1E answer, so I'm sure there's savings to be had here.

The magic square checking is borrowed from Kevin Cruijssen.

Explanation:

n      # Square input (implicit) (3 → 9)
 L     # Generate list from 1 to n^2 ([1,2,...,9])
  œ    # All permutations
   ʒ   # Filter by:
¹       # Recover n by pushing first input again
        # Check if magic square, borrowed from Kevin Cruijssen
 ô      # Split permutation into parts of size n
  D     # Duplicate
   ©    # Store in register (without popping)
    ø   # Zip rows to get columns
®       # Push from register
 Å     # Take main diagonal
®       # Push from register
 Å/     # Take anti diagonal
)       # Flatten stack into one list
 O      # Take sum (of each row/column/diagonal)
  Ë     # Check if all values are equal
     } # End filter (to get magic squares)
ʒ      # Filter magic squares by:
 ø      # Zip together magic square and input (implicit)
ε       # Map
 ¬       # Push the input again
  _      # Input equals 0 (to produce mask)
s        # Manage stack (swap mask and zipped args)
 Ë       # Partial equals potential match
  ~      # Bitwise OR to combine masks
    }   # End map
P      # Take product (effectively logical AND) to verify
       # that combined mask is all 1s
      # Implicit output

Brachylog, ⁴⁷ 44 bytes

{0∧|}ᵐ²{l⟦₅gj↔ʰc;?z∋₍ᵐġ,?;?ᵗc+ᵐ=&c≠≤ᵛ√~l?≜}

Try it online!

How it works

In classic Prolog style we replace zeros with uninitiated variables and based on the constraints let Brachylog figure out a solution. In Prolog you could just write [1,_,_] for unknown variables, in Brachylog you would have to write [1,A,B] and that seems too far away from the usual I/O restriction. So we use 0 for unknowns and convert them to uninitiated variables by:

{∧0|}ᵐ²

If a value is 0 try something else, otherwise use the value itself.

l⟦₅gj↔ʰc;?z∋₍ᵐ
l               length of array, N
 ⟦₅             0…N-1
   gj           [0…N-1],[0…N-1]
     ↔ʰc        0…N-1,N-1…0
        ;?z     [[0,first row], …, [N-1,last row],
                 [N-1,first row], …, [0,last row]]
           ∋₍ᵐġ [diagonal , diagonal /]

This feels a bit long just to get the two diagonals. Basically calculate the indices, zip them with the rows, and get the elements.

,?;?ᵗc

Append all rows and all transposed rows.

+ᵐ=

Sum every row. All sums (f.e. 15 in the 3x3 case) must be equal to each other. We don't have to calculate 15 explicitly, as this follows from the next constraint:

&c≠≤ᵛ√~l?
 c            the rows concatenated
  ≠           all elements are different
   ≤ᵛ         and are less-equal than X,
     √        and the root of X is
      ~l?     the length of the input
              which is implicitly the output

The numbers are distinct and between 1 and N^2.