Desmos, 94 92 90 bytes

f(n)=0^{0^{4^{round(u+.5log_4dd)-u}-x}}forx=L,y=L
d=3y-2^n2
u=(0^{0^d}+n)/2
L=[2^n...1]-.5

Try it on Desmos!

The matrix is returned in row-major order.

-2 bytes thanks to @Aiden Chow (log_4(dd)/2 → .5log_4dd)

-2 bytes now that list comprehensions don't require brackets on the RHS of assignments

Golfing note: 0^{0^{E-x} can be converted to \{E>x,0\}, but this requires moving E=round(u+.5log_4d)-u to its own line and doing \for instead of just for, so that is equal or +1 byte.

How it Works

The top third of matrix Bernard consists of right-aligned blocks of height \$2^{n-2},2^{n-4},\ldots,2^{\operatorname{mod}(n,2)}\$, and the bottom two-thirds consists of right-aligned blocks of height \$2^{n-1},2^{n-3},\ldots,2^{\operatorname{mod}(n+1,2)}\$. Each block has width equal to twice its height, and the top and bottom converge near y=m=2^n/3.

The main difficulty is for each y, find the height of the current block.

First we find the distance of y from the convergence y-value m. This is \$D=y-\frac{2^n}{3}\$

To figure out how many steps this is away, consider a situation where n is even, so the top third has step sizes \$1,4,16,\ldots\$. The partial sums of this geometric sequence are \$p(k)=\frac{4^k-1}{3}\$, e.g. \$p(2)=\frac{4^2-1}{3}= 5 = 1+4\$.

Then \$(x,y)\$ on the top third is in step size \$4^k\$ iff \$p(k)\leq D<p(k+1)\$. This can be re-arranged to \$4^k\leq 3D+1\leq 4^{k+1}\$, so \$k=\operatorname{floor}(\log_4(3D+1))\$

Ignoring golfing and the bottom side, this would give {x < floor(log_4(3D+1)): 1, 0} for the condition for filling in 1.

The final golfed code uses \$d=3y-2^{n+1}\approx 3D+1\$ (the y values are flipped, so this is really proportional to the distance of y from 2^n*2/3, explaining 2^{n+1} instead of 2^n).

Hence we have {x < floor(log_4(d)): 1, 0} as the current condition. See what it looks like on Desmos.

There's a few problems left: nothing is filled in for the bottom two-thirds, and the heights of the top third is twice as big as it should be for odd n.

Let's deal with the first problem first by just replacing log_4(d) with log_4(abs(d)) to handle negative d, or equivalently log_4(d^2)/2, which golfs to .5log_4dd. Then our condition becomes {x < floor(log_4(d)): 1, 0}. See on Desmos.

The issue now is that the heights of the top third are twice too big for odd n, and the heights of the bottom two-thirds are twice too big for even n. The solution comes down to handling off-by-half in the exponent of 4^{...} by adding and subtracting u=(0^{0^d}+n)/2 withing/outside the round. If d is positive (top third) and n is even, then u is an integer, so it doesn't affect the round. If d is negative (bottom two-thirds) and n is odd, then u is an integer, so it doesn't affect the round. Otherwise, n is a half-integer, so it causes round to behave like 1+floor, fixing the sizing issues.

Desmos, 80 bytes

f(n)=g(y,x)forx=L,y=L
g(y,x)=\{y>0:g(1-2y,2x-1),xx<1>-y,0\}
L=([1...2^n]-.5)/2^n

Try it on Desmos!

The matrix is returned in row-major order.

The recursive function g can be plotted implicitly to get the plot on the square (0,1)^2: Try it on Desmos!.

The key trick over my previous solution is to use recursion, which wasn't around at the time of that solution.

Python 3.8 (pre-release), 72 bytes

f=lambda n,z=0:[l:='0'*z+2**n*'1'][n:]or f(m:=n-1,z+2**m)[::-1]+2**m*[l]

Try it online!

BQN, 19 bytes

{(∾⟜×⟜⋆0¨∾˘⌽)⍟𝕩≍≍1}

Try it here!

-1 thanks to @att

Explanation

Input is x.

• (...)⍟𝕩≍≍1 execute x times on matrix [[1]]...
• 0¨⊸∾˘ row-wise prepend matrix of 0s
• 1¨⊸∾ column-wise prepend matrix of 1s
• ⌽ column-wise reverse

Example of how one iteration works:

0 1      # start
1 1

0 0 0 1  # 0¨⊸∾˘
0 0 1 1

1 1 1 1  # 1¨⊸∾
1 1 1 1
0 0 0 1
0 0 1 1

0 0 1 1  # ⌽
0 0 0 1
1 1 1 1
1 1 1 1
```

JavaScript (Node.js), 80 bytes

f=n=>n?f(n-1).reduce((v,e)=>[p=[...e.map(_=>0),...e],...v,p.map(_=>1)],[]):[[1]]

Try it online!

Python 3.8 (pre-release), 75 bytes

f=lambda n,g=[]:n and f(n-1,[0]*(p:=2**~-n)+g)[::-1]+p*[g+[1,1]*p]or[g+[1]]

Try it online!

Uiua, 19 bytes

⍥(⊂⊃∘=.⍉⊂⊃-⍉.⇌):¤¤1

Try it online!

How it works

⍥(⊂⊃∘=.⍉⊂⊃-⍉.⇌):¤¤1
                   ¤¤1 1x1 matrix containing a single 1 
                  :    flip the values on the stack, bringing n to the top of the stack
⍥                      repeat the following n times:
                ⇌      reverse the matrix's rows, flipping it vertically
               .       duplicate
           ⊃-⍉        apply subtraction from the duplicate and transpose in parallel
                       subtraction from the same matrix results in a 0 matrix
        ⊂             join the resulted matrices vertically
      ⍉               transpose back, reverting the previous transpose and making the join horizontal
      .               duplicate
   ⊃∘=                apply identity (keep, do nothing) and equality against the duplicate in parallel
                      equality against the same matrix results in a 1 matrix
  ⊂                   join

Python 3 with numpy, 77 76 bytes

from numpy import*
f=lambda n:n<eye(1)or kron(c_[:4:2,1:3],f(n-1)[::-1]+1)>1

Try it online!

• n<eye(1) produces [[True]] for n=0, the base case, and [[False]] for positive n, causing the expression after it to be used.
• c_[:4:2,1:3] produces \$ \begin{pmatrix}0&1\\2&2\end{pmatrix}\$, and is 1 byte shorter than the direct way of writing that matrix.
• The grid for n-1 is recursively obtained and flipped vertically, and then 1 is added, making it consist of 1s and 2s. kron then produces a copy of it in the top-right quadrant (>1 where the pre–+1 version had a 1), a doubled copy of it in each bottom quadrant (always >1), and zero in the top-left quadrant (never >1).

brev, 190 bytes

(define(b n)(define(e d v)(make-list(expt 2(- n d))v)) `(,@(e 2 `(,@(e 1 0),@(e 1 1))),@((over(append(e 1 0)(e 2 0)x))(b(- n 2))),@(e 1(e 0 1))))(define(b 1)'((0 1)(1 1)))(define(b 0)'((1)))

236 with whitespace added back in:

(define (b n)
  (define (e d v)
    (make-list
     (expt 2 (- n d)) v))
  `(,@(e 2 `(,@(e 1 0) ,@(e 1 1)))
    ,@((over (append (e 1 0) (e 2 0) x)) (b (- n 2)))
    ,@(e 1 (e 0 1))))

(define (b 1) '((0 1) (1 1)))
(define (b 0) '((1)))

Bonus! I like coming up with my own solution but here's my implementation of Bubbler's idea, at 135 bytes (displayed here with 15 unnecessary bytes of whitespace for readability)

(define (u n) ((fn `(,@((over `(,@((over 0) x) ,@x)) (reverse x)) ,@((over `(,@((over 1) x) ,@((over 1) x))) x))) (u (- n 1))))

(define (u 0) '((1)))

PARI/GP, 52 bytes

n->matrix(l=2^n,,i,j,j+2^#binary(bitxor(i-1,l\3))>l)

For input \$n\$, the number of \$1\$s in the \$i\$-th row (\$1\$-index) is \$2^{\#\operatorname{binary}(\operatorname{bitxor}(i-1,\lfloor2^n/3\rfloor))}\$, where \$\#\operatorname{binary}\$ means the number of the binary digits. This is because the binary representation of \$\lfloor2^n/3\rfloor\$ is \$(1010\dots)\$.

Attempt This Online!

Vyxal Ṁ, 20 bytes

0?(ṘnEn›ẋJ)E1vẋ0ÞṪṘ∩

Try it Online!

Charcoal, 26 bytes

↑1ＦＥＮＸ²ι«‖↓ＵＯ⊗ιι1Ｍ⊕ι↑»‖ＵＢ0

Try it online! Link is to verbose version of code. Explanation: Works by drawing the horizontal mirror image of Bernard, using the property that each Bernard is the vertical mirror image of the previous Bernard extended to the left with 0s and then extended down with 1s.

↑1

Output the lone 1 on its own row, leaving the cursor above the row.

ＦＥＮＸ²ι«

Loop n times, adding 2ⁱ rows each time, thus resulting in a total of 2ⁿ rows.

‖↓

Reflect vertically, leaving the cursor below Bernard.

ＵＯ⊗ιι1

Draw the next size rectangle of 1s.

Ｍ⊕ι↑

Move back to the top of the Bernard.

»‖

Reflect horizontally to complete the Bernard.

ＵＢ0

Fill the background with 0s.

Haskell, 111 bytes

data T a=C a|Q(T a)(T a)(T a)(T a)
f 0=C 1
f 1=Q(C 0)(C 1)(C 1)$C 1
f n=Q(C 0)(Q(C 1)(C 1)(C 0)$f$n-2)(C 1)$C 1

The matrix is represented as a quadtree. (I don't know if this output format is accepted).

Python, 121 104 bytes

lambda n:[[*map(int,f"{m:0{2**n}b}")]for m in g(n)]
g=lambda n:[1][n:]or g(n-1)[::-1]+2**~-n*[2**2**n-1]

Attempt This Online!

This is a much more interesting method, but is not well golfed because Python isn't great at converting binary numbers into a matrix (at least as far as I am aware).

Works by taking a start list of [1]. For each integer i up to n, reverse the list and add to it 2**~-n*[2**2**n-1]. This outputs the following pattern

[1, 3]
[3, 1, 15, 15]
[15, 15, 1, 3, 255, 255, 255, 255]
[255, 255, 255, 255, 3, 1, 15, 15, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535]

Which is then converted into binary, with each item in the list forming a row in the matrix.

-17 bytes thanks to @pxeger

Python, 119 117 bytes

f=lambda n:[[1]][n:]or-(-(p:=2**n)//4)*[(q:=p//2)*[0]+q*[1]]+[l+r for l,r in zip(p//4*[p//4*3*[0]],f(n-2))]+q*[[1]*p]

Attempt This Online!

Follows the algorithm as outlined in the question
-2 bytes thanks to pxeger

Jelly, 15 bytes

2Ṛ;»Ṁ²ƊƊ¡’Bz0ZU

A full program that accepts a positive integer from STDIN and prints a list representation of the Bernard.

Try it online! (The footer prints the result as a grid instead of in the raw list format.)

How?

Builds the Bernard in chunks consisting of the equal rows, starting with the row with the single 1. At each step we reverse what we have so far and append the next chunk. Each new chunk has as many rows as we've built so far but with double the number of 1s.

Under the hood, each row is first calculated as one more than its binary value. This starts at 2 and squares for each new chunk.

2Ṛ;»Ṁ²ƊƊ¡’Bz0ZU - Main Link: no arguments
2               - literal two
        ¡       - repeat STDIN times starting with X=2:
       Ɗ        -   last three links as a monad - f(X):
 Ṛ              -     reverse the current list
                        (2 gives its decimal digits reversed = [2])
      Ɗ         -     last three links as a monad - f(Z=that):
    Ṁ           -       maximum of Z
   »            -       vectorised maximum with Z -> an array of the same shape as
                                      Z, all populated with the maximum value in Z
     ²          -       square them
  ;             -     (reversed X) concatenate with (that)
         ’      - decrement all the values
          B     - convert them to binary lists
           z0   - transpose with filler zero  }
             Z  - transpose                   } prepends zeros
              U - reverse each

Python, 76 bytes

f=lambda n:(x:=1<<n>>1)and[x*[0]+j for j in f(n-1)[::-1]]+x*[x*[1,1]]or[[1]]

Attempt This Online!

Straight-forward recursion. Take the previous upside down, add a block of zeros on the left and a block of ones below.

Just seeing that @Bubbler used the same recursion a few hours before me.

05AB1E, 17 bytes

ÎF¼RNoF¾ª]oÅ10ζRø

Try it online or verify all test cases. (The footer is to pretty-print the matrix.)

Explanation:

Step 1: Create alternating lists of \$2^x\$ amount of 0-based values depending on the input:

Î          # Push 0 and the input
 F         # Loop `N` in the range [0, input):
  ¼        #  Increase `¾` by 1 (it's 0 by default)
   R       #  Reverse the current list (or 0 in the first iteration)
    No     #  Push 2 to the power `N`
      F    #  Inner loop that many times:
       ¾ª  #   Append `¾` to the list
           #   (`ª` will convert the 0 to [0] in the very first iteration)
 ]         # Close both loops

See just this step online for all test cases.

Step 2: Calculating \$2^v\$ for each of these values determines the amount of trailing 1s per row:

  o        # Map each value to 2 to the power that value
   Å1      # Then convert each inner value to a list of that many 1s

See the first two steps online for all test cases.

Step 3: Add leading 0s to each row, and output the resulting matrix:

      ζ    # Zip/transpose this list of rows; swapping rows/columns,
     0     # using a 0 as filler-character for unequal length rows
       R   # Reverse the list of lists
        ø  # Zip/transpose the matrix back
           # (after which it is output implicitly as result)

JavaScript (ES6), 95 bytes

A rather straightforward recursive construction.

n=>[...Array(1<<n)].map((_,y,a)=>a.map(g=(k=n,x,_,Y=y)=>k--?Y>>k|x>>k&g(k,x%=m=1<<k,_,~Y%m):1))

Try it online!

Commented

n =>                  // n = input
[...Array(1 << n)]    // build an array of 2 ** n entries
.map((_, y, a) =>     // for each item at position y in this array a[]:
  a.map(              //   iterate over a[] again using ...
  g = (               //   ... the recursive callback function g taking:
    k = n,            //     k = dimension, initialized to n
    x,                //     x = horizontal position
    _,                //     (ignored reference to the array)
    Y = y             //     Y = vertical position
  ) =>                //
    k-- ?             //     decrement k; if it was not 0:
      Y >> k |        //       insert a '1' if we are in the lower part (*)
      x >> k &        //       or we are in the right part and the result
      g(              //       of this recursive call is also 1:
        k,            //         pass the new value of k
        x %=          //         update x to x modulo m,
          m = 1 << k, //         where m = 2 ** k
        _,            //         (useless reference to the array)
        ~Y % m        //         update Y to (-Y - 1) modulo m
      )               //       end of recursive call
    :                 //     else:
      1               //       insert a '1' and stop the recursion
  )                   //   end of inner map()
)                     // end of outer map()

_{(*) The sign of Y is inverted at each iteration, so the 'lower' part is actually the upper part when Y is negative.}

Factor, 90 89 83 bytes

[ 2^ 3 dupn iota -rot '[ _ 3 /i bitxor bit-length 2^ 1 <array> _ 0 pad-head ] map ]

Try it online!

Port of @alephalpha's PARI/GP answer.

Python 3, 138 bytes

import re
def f(n):r=[int(bin(i)[2:])for i in range(2**n)];return[[1-bool(re.match('1(13)*(12|0)',str(x+y+y+10**n)))for x in r]for y in r]

Try it online!

K (ngn/k), 21 bytes

(|,/2\,'/4\10+)/[;+1]

Try it online!

Same algorithm, but uses a bit of base trick to generate the zero- and one-matrices at appropriate places.

10+    Map 0 to 10 and 1 to 11
4\     Base 4; 10 becomes (2 2) and 11 becomes (2 3)
2\     Base 2; 2 becomes (1 0) and 3 becomes (1 1)

K (ngn/k), 23 bytes

(,/(1|)\|,'/|^:\)/[;+1]

Try it online!

A straightforward implementation. Uses the recurrence pattern:

f(n+1) = [ zeros, f(n) reversed ]
         [ ones,  ones          ]

How it works

(,/(1|)\|,'/|^:\)/[;+1]    a function that takes n and
(               )/         applies the recurrence n times on
                    +1     the 1x1 matrix containing a single 1:
             ^:\    converge-scan using "test if each atom is null"; gives (mat; zeros)
         ,'/|       reverse the order and horizontally join the two matrices
        |           reverse vertically
   (1|)\      converge-scan using "max with 1"; gives (mat; ones)
 ,/           join vertically