I haven't used Matlab in a while, but something like this should work (for 100 points of p between 1 and 15, and for 50 points of q between 10 and 20):
p = linspace(1,15,100)
q = linspace(10,20,50)
for x = 1:length(p),
U[:,x] = ln(1./q) - cos(p[x])
end
plot3d(p,q,U)