There are codes on the FEX to compute the permanent of a matrix. One thing lacking is a tool with the capability to compute the permanent of a non-square matrix. This tool fills that gap. That was the ONLY reason I've posted this code, since there are already many tools on the FEX to compute a matrix permanent.
https://en.wikipedia.org/wiki/Permanent

If you want a reference on the algorithm employed, all I can suggest is:

https://en.wikipedia.org/wiki/Computing_the_permanent

The demos I've provided are pretty exhaustive. And since permanent works on symbolic arrays, I've offered afew tests for them too to convince you the computations are correct. But here are some simple examples of use:

permanent(magic(5))
ans =
53131650

permanent(rand(5,7) > 0.5)
ans =
181

A = sym('A',[2,3])
A =
[ A1_1, A1_2, A1_3]
[ A2_1, A2_2, A2_3]
permanent(A)
ans =
A1_1*A2_2 + A1_2*A2_1 + A1_1*A2_3 + A1_3*A2_1 + A1_2*A2_3 + A1_3*A2_2

For anyone who wants faster code, sorry. This is what it is. There are faster tools out there. permanent is pretty fast for matrices as large as 10x10, and is quite fast for smaller matrices. It does NOT use a recursive algorithm, but the fact is, it will need to compute the products of a lot of numbers for larger matrices.

tic,permanent(magic(10)),toc
ans =
4.73325614942977e+23
Elapsed time is 1.199238 seconds.

Computation for larger matrices can become VERY slow, because this tool will begin to use large blocks of memory. Recursive algorithms can also be slow for large problems.