or with MATLAB (the second result looks wrong to me):
syms x mu real
syms sigma real positive
syms k integer positive
f = exp(k*x)/(sqrt(2*sym(pi))*sigma)*exp(-(x-mu)^2/(2*sigma^2));
moment_k = int(f,x,-Inf,Inf)
double(subs(moment_k,[k mu sigma],[4 0.2 0.3]))
f1 = x^k/(sqrt(2*sym(pi))*sigma*x)*exp(-(log(x)-mu)^2/(2*sigma^2));
moment1_k = int(f,x,0,Inf)
double(subs(moment1_k,[k mu sigma],[4 0.2 0.3]))
