m-file:

tetha = pi/4;

lamb = -1;

h = 4;

tetha0 = 0;

syms x y l

n = [h.*((cos(tetha)).^2)./sin(tetha); h.*abs(cos(tetha)); 0];

ft = ((tetha - pi/2)./sin(tetha)).^4;

Rt = [cos(tetha) -sin(tetha); sin(tetha) cos(tetha)];

zt = [cos(tetha0) -sin(tetha0); sin(tetha0) cos(tetha0)];

lt = [x;y];

integrand = @(x,y)(ft.*h.*((abs(cos(tetha)).* (x.*cos(tetha)-y.*sin(tetha)))-((cos(tetha)).^2/sin(tetha)).*(x.*sin(tetha)+y.*cos(tetha))));

PhiHat = @(a,b)(dblquad(integrand,0,a,0,b));

ezsurfc(PhiHat,[0,5,0,5])