If you have the statistic toolbox, you'll want to use MVNCDF for a multivariate normal distribution. However, in case not everything you want to do is multivariate normal, let's fix your current approach. INTEGRAL, INTEGRAL2, and INTEGRAL3 assume that they can send in many U and K values at once and that you will return an array of corresponding p values. Yours doesn't work with non-scalar inputs. Assuming p = fun(U,K) works properly for all scalar U and K in the desired range of inputs, you can make it work with arrays using ARRAYFUN.
q = integral2(@(U,K)arrayfun(fun,U,K),1700,1900,30,120);
ARRAYFUN just takes the array inputs and loops over the entries to construct an array output.