H()

function out=H(W)

Input:

• W: Expansion coefficients for $N_s$ unconstrained wave functions, stored as an $\prod S_k\,\times\,N_s$ matrix

Global variables:

• gbl_Vdual: Dual potential coefficients stored as a $S_k\,\times\,1$ column vector.

Output:

• out: $H W$ (which represents ``LHS(W)'')

Using the expression in Eq. (11) and your Diagprod() operator above, produce a function of the above prototype in the file `H.m''