I would think the orthogonal array for 5-level factors would have to have a number of runs that is a multiple of 5, so it could not be 31.
In any case, there is no Statistics Toolbox function for computing general orthogonal arrays. There is a D-optimal function for generating designs for a specific combination of number of factors, number of levels, and number of runs. It usually gives good designs, and may reproduce the orthogonal array in some cases. Here's how you could generate a 31-run design for 4 5-level factors:
sortrows(rowexch(4,31,'linear','cat',1:4,'levels',5))
Sometimes when I run this, I get an orthogonal array with 25 runs:
sortrows(rowexch(4,25,'linear','cat',1:4,'levels',5,'tries',200))
