Plotting a Density Function Defined in Terms of an Integral

I want to plot the following function in 3D, but I first have to integrate w.r.t. y numerically since its not possible to get an analytical form.

Normally, I woulld just use

x_1 = linspace(-a,a,300);

x_2 = linspace(-a,a,300);

[X, Y] = meshgrid(x_1,x_1);

but the integral makes things a bit difficult. How can I work around this?