How to efficiently calculate matrix operations when the elements are defined as functions?

I'm trying to calculate the one-body density matrix of a system of particles as described in this article, given in equations 14 and 15.

Here, is my code but it seems very slow. Is there any efficient way to calculate it in Mathematica?

\[Omega]z = 2 \[Pi] 14.3;
hbar = 1.05457*10^-34; (* Planck's constant *)
kb = 1.38065*10^-23; (* Boltzman's constant *)
mCs = 132.9054*1.6605*10^-27; (* mass Cs-Atom *)
lho = Sqrt[hbar/(mCs *\[Omega]z)];
Lmax = 10;
Np = 5;
\[Mu] = 2.3*10^-32;

Energy[n_] := (2 n + 1) *(hbar/2)* \[Omega]z;

f[i_] := 1/(Exp[(Energy[i] - \[Mu])/(kb*T)] + 1);

\[Psi][z_][n_] :=
N[Exp[-0.5 (z/lho)^2]*
HermiteH[n, z/lho]*(1/(N[\[Pi]]^(1/4) Sqrt[2^n n! lho]))] ;
\[Psi]bar[z_][n_] :=
N[Exp[-0.5 (z/lho)^2]*
HermiteH[n, z/lho]*(1/(N[\[Pi]]^(1/4) Sqrt[2^n n! lho])) /.
I -> -I ];

Pjk[z_][j_, k_] :=
Pjk[z][j, k] =
Exp[-z^2]*(HermiteH[j + 1, z]*HermiteH[k, z] -
HermiteH[j, z]*HermiteH[k + 1, z])*(1/(
2*(j - k)*Sqrt[2^(j + k) \[Pi] j! k!]));
Mj[z_, 0] = 0;
Mj[z_, j_] :=
Mj[z, j] =
Mj[z, j - 1] +
Exp[-z^2]*HermiteH[j, z]*
HermiteH[j + 1, z]*(1/(2^(j + 1)*((j + 1)!))) // N;
P[x_, y_][j_, k_] :=
P[x, y][j, k] =
If[j == k,
1 - 2*Sign[y - x]*f[j - 1]* (Mj[y, j - 1] - Mj[x, j - 1]), -2*
Sign[y - x]*
Sqrt[f[j - 1]*f[k - 1]]*(
Pjk[y][j - 1, k - 1] - Pjk[x][j - 1, k - 1])];
pmat[x__, y__] := Array[P[x, y], {Lmax, Lmax}] // N;
(*qmat[x_,y_]:=ConjugateTranspose[Inverse[pmat[x,y]]]*Det[pmat[x,y]];*)

qmat[x_, y_] :=
N[ConjugateTranspose[LinearSolve[pmat[x, y], IdentityMatrix[Lmax]]]*
Det[pmat[x, y]]];
allpsi[x_] :=
allpsi[x] = Table[\[Psi][x][i]*Sqrt[f[i]], {i, 0, Lmax - 1}] // N;
allpsibar[y_] :=
allpsibar[y] =
Table[\[Psi]bar[y][i]*Sqrt[f[i]], {i, 0, Lmax - 1}] // N;
densitymat[x_, y_] :=
Parallelize[allpsi[x] . qmat[x, y] . allpsibar[y]];
Plot[densitymat[x, 0], {x, 0, 5 lho}, PlotRange -> All]


For example, for (Lmax=10, Np=5), it takes around 5s, for(Lmax=20,Np=10), it's around 80s, and for (Lmax=50, Np=10), it's more than 600s. Is there any way to perform these matrix operations more efficiently so that we can push it to higher values of (Lmax,Np)? Finally, I would like to do a Fast Fourier transform of densitymat.