What do you do in Mathematica? What is your goal in the end?
The transverse-field Ising model you solve above could be implemented in Julia rather easily as well. I study disordered spin systems and have some libraries for constructing spin models and stuff. You could write an equivalent code to your OP like:
using SpinModels # not registered - see https://github.com/abraemer/SpinModels.jl
function psi_t(N, h, t)
# TFIM
H = NN(PBC(Chain(N))) * ZZ(0.5) + h*X(ones(N))
U = cis(Hermitian(Matrix(-t*H)))
state = zeros(2^N)
state[1] = 1.0
return U*state
end
Of course this is not efficient if you want different times and different states. So what is it you want to achieve in the end
Also if you wanted to exploit some symmetries, you can use SpinSymmetries.jl to construct projectors onto the differen sector. This is usually still worth it even if you state lives in all sectors (as it does in this example)