Desperate for help. Calculate VIX option price with Grünblicher & Longstaff (1996)

Please someone help. I have build the below function to price a Call option on the VIX. I have directly replicated the call price function from Grünblicher & Longstaff (1996), but I am still not able to reproduce the values they achieve in Fig. 1 in their article.

____________________________________________

function F=hestonsv(theta,K,T,r,VIX)

bet=theta(1);   lrv=theta(2);   sig=theta(3);

gam=(4.*bet)./((sig^2).*(1-exp(-bet.*T)));

vega=(4.*lrv)./(sig^2);

lambda=gam.*exp(-bet.*T).*VIX;

F=exp(-r.*T).*exp(-bet.*T).*VIX.*ncx2cdf(gam.*K,(vega+4),lambda)+exp(-r*T).*(lrv./bet).*(1-exp(-bet.*T)).*ncx2cdf(gam.*K,(vega+2),lambda)-exp(-r.*T).*K.*ncx2cdf(gam.*K,vega,lambda);
___________________________________________________

I am using the exast same input parameters as in their article. r=0.05, lrv=0.6, bet=4, sig^2=0.133 and K=0.15. No matter how i twist and turn it i cannot reproduce their result. Someone please help, what am i doing wrong?