% computes the hessian of F using the differential approximation to |x|
function H=hessianF(A,At,b,p,lambda,sigma,x)

    m = size(A,1); n = size(A,2);
    erf_vec = zeros(n,1);
    alpha = sqrt(2/pi)*(sigma);
	hes_mat1 = 2*At*A;
	hes_mat2 = eye(n,n);

    for j=1:n
        xj = x(j);
		erf_val = erf(xj/(sqrt(2)*sigma));
		phi_sig_xj = xj*erf_val + alpha*exp(-xj^2/(2*sigma^2));
		hes_mat2(j,j) = (p-1)*(phi_sig_xj)^(p-2) * erf_val^2 + (phi_sig_xj^(p-1))*sqrt(2/(pi*sigma^2)) * exp(-xj^2/(2*sigma^2));
    end
    H = hes_mat1 + lambda*p*hes_mat2;

end
