plotting heat diffusion in a sphere

I am trying to plot heat diffusion in a sphere but I am having no luck. I am trying to adapt my code for a square.

My code for a square is:

% Initialize square grid
[x, y] = meshgrid(0:0.005:1, 0:0.005:1);

nmax = 110;

T = zeros(length(x), length(y));

for n = 1:nmax
    m = 2*n - 1;
    T = T + 400/pi*sinh(m*pi*y)/(m*sinh(m*pi)).*sin(m*pi*x);
end

surface(x, y, T)
shading('interp')
colorbar
colormap jet

For the sphere, I have written (which doesn't work):

[x, y] = meshgrid(-1:.01:1, -1:.01:1);

[theta, r] = cart2pol(x, y); 
N = 3;

p = zeros(length(N));
T = zeros(length(x), length(y));

fun = @(x) legendre(0, x);

for l = 0:N
    p = p + 50*(2*l + 1).*integral(fun, 0, 1);
end

for n = 1:N
    T = T + p.*r.^n.*legendre(0, cos(theta));
end

surface(theta, r, T) shading('interp') colorbar colormap jet

The Legendre integral isn't returning the correct coefficients. The first 3 are 50, 75, 0. The code is returning 50, 200, 450.

What is going wrong?