$\newcommand{\ket}[1]{\lvert#1\rangle}$
I don't know why the methods are more suited for humans than for computer, and you can easily find a MatLab implementation by googling, for example on Toby Cubitt's webpage. (I haven't checked the algorithm, but I trust Toby to get such a thing correctly.)
I'll try below to give you an explanation “suited for computers”. Let $\rho_{AB}$ be the density matrix you're interested in, with the $d_A d_B$ base vectors in the following order :
$$\ket{1,1},\ket{1,2},\dots,\ket{1,d_B},\ket{2,1},\dots,\ket{d_A,d_B}$$
To obtain $\rho_B=\text{Tr}_A(\rho_{AB})$, you just have to add the $d_A$ diagonal submatrices of $\rho_AB$. Something like $\rho_B=\sum_{j=0}^{d_A-1}\rho_{AB}[jd_B+1\dots jd_B+d_A;jd_B+1\dots jd_B+d_A]$.
To get $\rho_A$, you can reorder the indices to go back to the previous algorithm. Or you directly take the relevant submatrices, but they are somehow scattered. $\rho_A[k;l]=\sum_{i}^{d_B}\rho_{AB}[(k-1)d_B+i;(l-1)d_B+i]$