How to make a surface in polar coordinates using polar3d?

This takes too long to run here (it took 336.315050 seconds — 00:05:36.315049 — just now on MATLAB Online) however it plots the surface on a ‘synthetic’ polar coordinate axis in 3D. It would be relatively straightforward to change the Z axis coordinates (or ellimiinate them entirely). The original problem was that the code for the polar axes was not complete. It now works, as it worked in my earlier code.

Try This:

s = 3;
n = 1;
p = 0.1; % this is time
tv = 0:10:360; %3; % this is angle theta for polar rotation
r = 1;
a = 0:10:360;
a = deg2rad(a);
tv = deg2rad(tv);
tic
for k = 1:numel(tv)
    t = tv(k)
    b = sqrt(2*n*p);
    L = sqrt((4*p+r^2)/3);
    fun = @(k,u,c,a) ((k.^2).*exp(-1.5*k.^2)).*((u.^2).*(1-u.^2).*exp(-(b.*u).^2).*(cos(s.*k.*u.*cos(a)/L))).*(((cos(c)).^2).*(cos(k.*sqrt(1-u.^2).*(s.*sin(a).*cos(t).*cos(c)+s.*sin(a).*sin(t).*sin(c))/(L)))); 
    f3 = arrayfun(@(a)integral3(@(k,u,c)fun(k,u,c,a),0,Inf,-1,1,0,2*pi),a);
    B = ((6*sqrt(6)*(b^3))/(erf(b)*(pi^2)))*(1-(3/(2*b^2))*(1-((2*b*exp(-b^2))/(erf(b)*sqrt(pi)))))^(-1);
    R = B*f3;      
    ta = t*ones(size(tv));
    [x,y,z] = pol2cart(a, R, ta);
    X(k,:) = x;
    Y(k,:) = y;
    Z(k,:) = z;
    toc                            % Total Time To This Point
    LIT = toc/k                    % Mean Loop Iteration Time                       
end
toc
figure (2)
surfc(X, Y, Z)
colormap(turbo)
axis('equal')
axis('off')
hold on
a = deg2rad(0:2:360);               % Create  & Plot Polar Cylindrical Coordinates
r = 1;
xc = [0:0.25*r:r].'*cos(a);
yc = [0:0.25*r:r].'*sin(a);
xr = [0;r]*cos(0:pi/4:2*pi);
yr = [0;r]*sin(0:pi/4:2*pi);
plot3(xc.', yc.', zeros(size(xc)), '-k')
plot3(xr, yr, zeros(size(xr)), '-k')
zt = (0:1:max(tv)).';
zt1 = ones(size(zt));
surf(zt1*xc(end,:), zt1*yc(end,:), zt*ones(size(a)), 'FaceAlpha',0, 'MeshStyle','row')
hold off
text(xr(2,1:end-1)*1.1,yr(2,1:end-1)*1.1, zeros(1,size(xr,2)-1), compose('%3d°',(0:45:315)), 'Horiz','center','Vert','middle')
text(ones(size(zt))*xc(end,end-10), ones(size(zt))*yc(end,end-10), zt, compose('%.1f',zt))
toc
TTmin = toc/60