Filling out an $n \times n$ square grid with  $0$s and $1$s 
How many ways are there of filling an $n×n$ square grid with $0$s and $1$s
  if you are allowed at most two $1$s in each row and two $1$s in each
  column?

I need some ideas for solving this problem.
PS:This problem is from the book Princeton Companion Mathematics.
 A: Due to popular demand, here's an explanation of the recurrence used in the code I posted at OEIS.
Let $a(k,l,m,n)$ be the number of binary $m\times n$ matrices with at most two $1$s per row, $k$ column sums of $0$, $l$ column sums of $1$ and the remaining $n-k-l$ column sums of $2$. In the following, I'll call a column with sum $n$ an $n$-column for short. Since the numbers we're looking for have no restrictions on $k$ and $l$, they're given by 
$$a(n)=\sum_{k=0}^n\sum_{l=0}^{n-k}a(k,l,n,n)\;.$$
There is exactly one $0\times n$ matrix, and it has $n$ column sums of $0$, so $a(n,0,0,n)=1$ and $a(k,l,0,n)=0$ otherwise. To find a recurrence relation for $a(k,l,m,n)$, we can go through all admissible ways of adding a new row to an $(m-1)\times n$ matrix.
We can add a row without $1$s; that doesn't change the column sums and gives a contribution $a(k,l,m-1,n)$.
We can add a row with one $1$, and the $1$ can be in a $0$-column or a $1$-column. If it's in a $0$-column, it reduces $k$ and increases $l$, and there are $k+1$ columns to choose from, for a contribution of $a(k+1,l-1,m-1,n)(k+1)$. If it's in a $1$-column, it leaves $k$ unchanged and reduces $l$, and there are $l+1$ columns to choose from, for a contribution of $a(k,l+1,m-1,n)(l+1)$.
Or we can add a row with two $1$s, and the columns they're in can be either two $0$-columns, one $0$-column and one $1$-column, or two $1$-columns. If they're two $0$-columns, $k$ is reduced by $2$ and $l$ is increased by $2$, and there are $(k+2)(k+1)/2$ pairs of columns to choose from, for a contribution of $a(k+2,l-2,m-1,n)(k+2)(k+1)/2$. If they're one $0$-column and one $1$-column, $k$ is reduced by $1$ and $l$ is unchanged, and there are $(k+1)l$ pairs of columns to choose from, for a contribution of $a(k+1,l,m-1,n)(k+1)l$. If they're two $1$-columns, $l$ is reduced by $2$ and $k$ is unchanged, and there are $(l+2)(l+1)/2$ columns to choose from, for a contribution of $a(k,l+2,m-1,n)(l+2)(l+1)/2$.
Summing it all up, we have
$$
\begin{eqnarray}
a(k,l,m,n)
=&&
a(k,l,m-1,n)
\\
&+&a(k+1,l-1,m-1,n)(k+1)+a(k,l+1,m-1,n)(l+1)
\\
&+&
a(k+2,l-2,m-1,n)(k+2)(k+1)/2
\\
&+&
a(k+1,l,m-1,n)(k+1)l
\\
&+&
a(k,l+2,m-1,n)(l+2)(l+1)/2\;,
\end{eqnarray}
$$
where terms with indices out of range are taken to be $0$.
