AWK, 67 bytes

{for(i=int(sqrt($1)+1);$1%--i;);for(x=$1/i;i--;)printf"%0*d\n",x,0}

Attempt This Online!

{for(i=int(sqrt($1)  # start nearest to square
+1);                 # better square shape
$1%--i;);            # variation on primality test
for(x=$1/i           # get other side length
;i--;)               # count back down 
printf"%0*d\n",x,0}  # print string of 0s x length

SAKO, 127 126 bytes

CZYTAJ:N
*)GDYMOD(ENT(N),ENT(A))=0:1,INACZEJ2
2)POWTORZ:A=ENT(PWK(N))(1)N
*1)SPACJIENT(N/A)
LINIA
POWTORZ:I=1(1)A
STOP1
KONIEC

Prints a A by N/A grid of spaces.

CASIO BASIC (CASIO fx-9750GIII), 84 bytes

?→A~D
For Int √A→J To .5A
A÷J=A Int÷ J And J-A˼J<D-A˼D
DNot Ans+JAns→D
Next
For 1→Q To D
For 1→R To A÷D
PxlOn R,Q
Next
Next

I apologize for the wait.

it draws the almost-square in graph mode

C (gcc), 60 bytes

j,k;f(n){for(j=k=n;j--;)j*j>n|n%k?k--:printf("%0*d\n",n/k);}

Try it online!

Standard math library (i.e. sqrt) not used. Instead, candidate integers are iterated over in descending order, starting from the input number. Integers are tested for whether they are a divisor of the input number and whether they are near or below the square root. Once such a number is found, the rectangle is printed.

Similar approach is used above for JS and Python. This one was derived from the other C solution posted by jdt.

By some stroke of luck, it appears that there is a zero in the third argument on TIO, which saves two bytes in the call to printf. Consequently the function may randomly fail for different callers.

Ternary conditional results in a 1-4 byte reduction from other golfy approaches. I suppose that I could test whether the iterator is at the square root exactly, but where's the fun in that.

05AB1E, 11 10 9 bytes

ÑÅsD‚˜`ии

Outputs as a matrix of either the input (for squares) or one of the two side-lengths (for non-squares).

Try it online or verify all test cases.

Explanation:

          #  Example input 1: 12
          #  Example input 2: 16
Ñ         # Get the divisors of the (implicit) input-integer
          #  STACK 1: [1,2,3,4,6,12]
          #  STACK 2: [1,2,4,8,16]
 Ås       # Get the middle value (for squares) or middle pair (for non-squares)
          #  STACK 1: [3,4]
          #  STACK 2: 4
   D‚     # Pair it with a copy of itself
          #  STACK 1: [[3,4],[3,4]]
          #  STACK 2: [4,4]
     ˜    # Flatten, so we now have a pair or quartet
          #  STACK 1: [3,4,3,4]
          #  STACK 2: [4,4]
      `   # Pop and push both or all four values to the stack
          #  STACK 1: 3,4,3,4
          #  STACK 2: 4,4
       и  # Repeat the second-last item the last item amount of times as list
          #  STACK 1: 3,4,[3,3,3,3]
          #  STACK 2: [4,4,4,4]
        и # Map each inner value to a list of that many time the top value (for
          # non-squares), or the implicit input (for squares)
          #  STACK 1: [[4,4,4],[4,4,4],[4,4,4],[4,4,4]]
          #  STACK 2: [[16,16,16,16],[16,16,16,16],[16,16,16,16],[16,16,16,16]]
          # (after which the matrix is output implicitly as result)

NOTE: I can't remove the ¸˜, since the ` would incorrectly push the loose digits of multi-digit squares (e.g. input=\$100\$ with middle divisor \$10\$ would push 1,0 instead of 10,10 without the ¸˜).

APL(Dyalog Unicode), 29 24 18 bytes ^SBCS

Find the middle factor pair, then reshape the number into the rectangle.

18 bytes

0⍴⍨∘⊃⊢⍸⍤=⍳∘.(××≥)⍳
0⍴⍨∘⊃              ⍝ 0 in the shape of the first of
     ⊢⍸⍤=⍳∘.(××≥)⍳ ⍝ the factor pairs that its left is not less than its right.
                   ⍝ More precisely,
     ⊢⍸⍤=          ⍝ indices of N in
         ⍳∘. ×   ⍳ ⍝   the mulultiplication table of [1, N]
            ( ×≥)  ⍝     with its upper-right triangle being set to 0s

Try it on APLgolf!

24 bytes

{⍵⍴⍨i⊃⍨⌈2÷⍨⍴i←⍸⍵=∘.×⍨⍳⍵}
                 ∘.×⍨⍳⍵  ⍝ Multiplication table of [1,⍵]
            i←⍸⍵=        ⍝ Factor pairs with product being ⍵
    i⊃⍨⌈2÷⍨⍴             ⍝ The pair in the middle
 ⍵⍴⍨                     ⍝ ⍵ in the shape of the pair   

Try it on APLgolf!

29 bytes

(⊣⍴⍨⊢,÷)∘((⍸⊃⍨2⌈⍤÷⍨+⌿)0=⍳|⊢)⍨

It's the tacit version of {⍵⍴⍨j,⍵÷j←(⍸i)⊃⍨2⌈⍤÷⍨+⌿i←0=⍵|⍨⍳⍵}. The following snippet explains what it does more clearly.

{
    i ← 0=(⍳⍵)|⍵      ⍝ Bit mask for factor of N in [1,N]
    j ← (⌈(+⌿i)÷2)⊃⍸i ⍝ The factor in the middle
                      ⍝ (Take the larger one for odd N)
    k ← j,⍵÷j         ⍝ The desired shape
    k ⍴ ⍵             ⍝ Reshape of N in the desired shape
}

Try it on APLgolf!

25 byets with sqare root

Another one using square root in 25 bytes, but I like the other ones without it more.

{⍵(⊣⍴⍨⊢,÷)0⊥⍸0=⍵|⍨⍳⌈⍵*.5}

Try it on APLgolf!

Vyxal, 7 bytes

K:Y∆ṁÞm

Try it Online!

Explanation:

K gives the increasing list of all of the factors of the input.

:Y interleaves this list with itself, giving two of each factor.

∆ṁ gives the median of the list. The list is always of even length because we doubled the number of elements previously.

In most cases we wouldn't need to do the interleaving step, as the two middle elements would be the factor pair. In the case of the input already being a perfect square, however, that number is the middle of the list but it's only present once. By having two of each factor, we fix this case, and the normal case is unaffected because the two middle elements are still the factor pair.

Finally, Þm turns the list of dimensions into an array with that shape filled with zeroes.

C (gcc), 66 bytes

j;f(n){for(j=sqrt(n);n%j;j--);for(n/=j;j--;)printf("%0*d\n",n,0);}

Try it online!

Google Sheets, 78 bytes

=let(i,sequence(A1^0.5),y,max(filter(i,0=mod(A1,i))),rept(rept(0,A1/y)&"
",y))

Put \$N\$ in cell A1 and the formula in cell B1. The formula uses a hard-coded newline in quotes — you can enter that with Control + Enter while editing the formula.

The result is a text string with newlines. To get a 2D array instead, use this (85 bytes):

=let(i,sequence(A1^0.5),y,max(filter(i,0=mod(A1,i))),makearray(y,A1/y,lambda(r,c,y)))

This formula uses the number of rows as the value to show in the grid.

Ungolfed:

=let( 
  candidates, sequence(A1 ^ 0.5), 
  height, max(filter(candidates, not(mod(A1, candidates)))), 
  width, A1 / height, 
  makearray(height, width, lambda(row, column, height)) 
)

(-5 bytes thanks to jdt)

JavaScript (Node.js), 52 bytes

f=(n,q=r=`
`)=>q*q<n|n%q?f(n,++q,r+=0):r.repeat(n/q)

Try it online!

From Arnauld's

Python 3, 50 bytes

f=lambda n,q=1:q*[n//q*[0]][n%q:q*q>=n]or f(n,q+1)

Try it online!

-3 bytes from xnor and Albert.Lang

Same idea

Perl 5 -pa, 45 bytes

$;=0|sqrt;$;++while$_%$;;$_=(1x$;.$/)x($_/$;)

Try it online!

JavaScript (ES6), 50 bytes

Returns a single string made of 0's and newlines.

f=(n,r=q=`0
`)=>++q*q<n|n%q?f(n,0+r):r.repeat(n/q)

Try it online!

Commented

f = (             // f is a recursive function taking:
  n,              //   n = input
  r =             //   r = string holding the the almost-square row,
                  //       initially set to "0" + newline
  q = `0\n`       //   q = divisor, initially zero'ish
) =>              //
++q * q < n |     // increment q; if q² is less than n
n % q ?           // or q is not a divisor of n:
  f(              //   do a recursive call:
    n,            //     pass n unchanged
    0 + r         //     insert a "0" at the beginning of r
  )               //   end of recursive call
:                 // else:
  r.repeat(n / q) //   return r repeated n/q times

Charcoal, 16 bytes

Ｎθ⪪×0θ⊕⌈Φ₂θ¬﹪θ⊕ι

Try it online! Link is to verbose version of code. Only works up to 2¹⁰⁶ due to floating-point inaccuracy but but would be stupidly slow well before then anyway. Outputs a rectangle of 0s where the width is no greater than the height. Explanation:

Ｎθ                  Input `N` as a number
    0               Literal string `0`
   ×                Repeated
     θ              `N` times
  ⪪                 Split into substrings of length
          θ         Input `N`
         ₂          Square root
        Φ           Filtered over implicit range where
             θ      Input `N`
           ¬﹪       Is divisible by
               ι    Current value
              ⊕     Incremented
       ⌈            Take the maximum
      ⊕             Incremented
                    Implicitly print

Retina 0.8.2, 46 bytes

.+
$*
^((.)+?)((?<-2>\1))*$
$1$#3$*0
+1`0
¶$%`

Try it online! Link includes test cases. Outputs a rectangle of 1s (this can be changed without affecting the byte count) where the height is no more than 1 greater than the width. Explanation:

.+
$*

Convert to unary.

^((.)+?)((?<-2>\1))*$
$1$#3$*0

Factorise the input into two almost-square-roots. $#3 is actually 1 less than the second factor, which is why it might equal $.1 and therefore the second factor would be 1 greater, although they are still both almost square roots in this scenario. The first factor is represented in unary as a string of 1s while the second factor can use any other printable ASCII character as it will be replaced with newlines below. Alternatively, the first character could have been any ASCII character other than 1 in which case the second factor would need to be represented by 1s to keep the same byte count.

+1`0
¶$%`

Repeat the first factor the second factor of times.

38 bytes in Retina 1:

.+
*
^((.)+?)((?<-2>\1))*$
$1$#3*$(¶$1

Try it online! Link includes test cases. Outputs a rectangle of _s where the height is no more than 1 greater than the width. Explanation: As above but $( is used to generate the output rectangle directly.

35 bytes in Retina 1 with a trailing newline:

.+
*
(?<-2>(^(_)+?|\1))+$
$#1*$($1¶

Try it online! Outputs a rectangle of _s where the height is not greater than the width. Explanation: Based on my answer to Closest to a square which factorises the input into $#1 and $.1 but outputs the rectangle instead of its dimensions.

TinyAPL, 18 bytes

÷⇾⍪⟜{⊇⍒⍵|⍨⍳⧺⌊√⍵}⊸⍴

Try it

The method this uses comes from nyxbird's Uiua solution.

I did some experimentation to find the shortest way to arrange all the pieces here, and I ended up with having the middle section written as an explicit function because writing it tacitly required a lot of ⍤ Atops and ⍛ Befores:

÷⇾⍪⟜(⍳⍤⧺⍤⌊⍤√⊸|⍛⍒⍛⊇)⊸⍴

Uiua, 18 14 bytes

^{2 away... birdie! -2 bytes thanks to noodle person}

↯⊟◡÷⊣⍖⊸◿⇡+1⌊√.

Try it!

↯⊟◡÷⊣⍖⊸◿⇡+1⌊√.
          ⇡+1⌊√   # get the range from 0-flr(sqrt(n))+1
         ◿     . # mod n by each
    ⊣⍖          # find the index (i) of the last 0
↯⊟◡÷  ⊸         # reshape n into a i x (n/i) array

⊣ last ⍖ fall gets the index of the smallest element, breaking ties toward the end. Here, 0 is always the smallest because at least one number in the range divides n.

Jelly,  8  7 bytes

sⱮ½¬EƇṪ

A monadic Link that accepts a positive integer, \$N\$, and yields a rectangular list of lists of zeros as near to square as possible (more rows than columns if non-square \$N\$). (Won't work for huge inputs due to floating point inaccuracies creeping in.)

Try it online!

How?

sⱮ½¬EƇṪ - Link: positive integer, N
  ½     - square root {N}
 Ɱ      - map across {x in [1 .. floor(root(N))]} with:
s       -   split {[1..N]} into chunks of length {x} (trailing shorter if need be)
   ¬    - logical NOT (converts all values to zeros)
     Ƈ  - keep those for which:
    E   -   all equal?
      Ṫ - tail

8 bytes

(With no floating point inaccuracy)

Either change ½ to ƽ in the above or...

sÆDŒHḢṪ¬

Try it online!

How?

sÆDŒHḢṪ¬ - Link: positive integer, N
 ÆD      - divisors of {N}
s        - {[1..N]} split into chunks of {those}
   ŒH    - split into two (with first half one longer when N is square)
     Ḣ   - head of {that} -> all ways to make a rectangle with rows <= root(N)
      Ṫ  - tail of {that} -> the one with the most rows
       ¬ - logical NOT (convert all the values to zeros)

Python 2, 74 bytes

def f(b):a=[x+1for x in range(b)if-1<b%~x<b>x*x][-1];print(b/a*'0'+'\n')*a

Try it online!

Japt, 7 bytes

Uses q, outputs an array of lines

ç¬òUâ Ì

Try it

ç¬òUâ Ì     :Implicit input of integer U
ç           :Repeat U times
 ¬          :  "q"
  ò         :Partitions of length
   Uâ       :  (Paired) Divisors of U
      Ì     :  Last element

Japt -h, 7 bytes

Uses l, outputs an array of lines

â £îÊòX

Try it

â £îÊòX     :Implicit input of integer U
â           :Divisors
  £         :Map each X
   î        :  Mold to length U
    Ê       :    "l"
     òX     :  Partitions of length X
            :Implicit output of last element

Wolfram Language (Mathematica), 47 46 40 bytes

Table[,If[i∣#,d=i]~Do~{i,√#};d,#/d]&

Try it online!

Returns an array of Nulls.

                  ~Do~{i,√#}        for i=1... i*i<=N:
       If[i∣#,d=i]                    i divides N: update shape
Table[,                      d,#/d] Nulls of that shape

R, 43 bytes

\(n,i=1:n^.5)matrix(0,k<-max(i[!n%%i]),n/k)

Attempt This Online!

Returns a 2D matrix of 0s.

R, 48 bytes

\(n,i=1:n^.5)write(rep(0,n),1,max(i[!n%%i]),,"")

Attempt This Online!

First version when I didn't read that printing output is optional not required.

APL+WIN, 36 34 bytes

Prompts for n

((n÷v),v←↑⌽(0=m|n)/m←⍳⌊(n←⎕)*.5)⍴0

Try it online! Thanks to Dyalog Classic