You can evaluate molecular orbitals on a grid by first evaluating the AOs on the grid with

 ao = mol.eval_gto('GTOval', coords)

and then contracting with the orbital coefficients $$ \psi_{i\sigma} = \sum_\alpha C^\sigma_{i\alpha} \chi_\alpha $$

I don't remember off-hand which way the indices are in PySCF, but this should be easy to figure out. The code should look something like

mo = lib.einsum('gi,xij->xgj', ao, mf.mo_coeff)

You can evaluate molecular orbitals on a grid by first evaluating the AOs on the grid with

 ao = mol.eval_gto('GTOval', coords)

and then contracting with the orbital coefficients $$ \psi_{i\sigma} = \sum_\alpha C^\sigma_{i\alpha} \chi_\alpha $$

I don't remember off-hand which way the indices are in PySCF, but this should be easy to figure out. The code should look something like

mo = lib.einsum('gi,xij->xgj', ao, mf.mo_coeff)

addendum regarding edit above:

ao vals shape:  (79368, 46)
mo vals shape:  (2, 79368, 46)
mo coeff shape:  (2, 46, 46)

This means you have 79368 grid points, 46 basis functions, and 46 spatial molecular orbitals (in general you may have fewer molecular orbitals than atomic orbitals as the atomic orbital basis set may have linear dependencies). For unrestricted wave functions, PySCF uses arrays with additional indices for the spin; I believe [0] is spin-up and [1] is spin-down.