Real spherical harmonics for signal approximation with high degree yield erroneous result

Question

I`m not sure this question fit to this area. If not, please let me know.

For reconstructing a signal $f(\theta,\phi)$ in terms of "Real" spherical harmonic basis $Y_{l,m}(\theta,\phi)$, i.e.,

$f(\theta,\phi) = {\sum}f_{l,m}Y_{l,m}(\theta,\phi)$,

the reconstructed signal gets weird with increasing $l$(degree).

I have tested with $f(\theta,\phi)=\max(0,\cos(\theta))$. The original signal on sampled sphere looks like this.

With $l=4$, the reconstructed signal seems reasonable to me.

However, when I grow $l$ to 10 it gets worse than smaller degree.

Also, the mean squared error with original signal is getting worse with growing number of degree $l$. Is this problem comes from limitation in computational bound? or any explanation for this phenomenon?

[Edit] : Thank you lisyarus, I attach a Matlab code to reproduce above result. I referred code example from "Green, Robin. "Spherical harmonic lighting: The gritty details." Archives of the Game Developers Conference. Vol. 56. 2003.",

test.m

clear all; close all;
[sph, vec, coef] = SH_setup_spherical_samples(10000, 4);
figure, scatter3(abs(coef(:,2)).*vec(:,1), abs(coef(:,2)).*vec(:,2), abs(coef(:,2)).*vec(:,3), [], coef(:,2));
colormap jet;
fn = @(theta, phi) max(0, 5*cos(theta)-4) + max(0, -4*sin(theta-pi).*cos(phi-2.5)-3);
results1 = SH_project_polar_function(fn, sph, coef);
fval = fn(sph(:,1), sph(:,2));
figure, scatter3(fval.*vec(:,1), fval.*vec(:,2), fval.*vec(:,3), [], fval);
colormap jet;

SH_setup_spherical_samples

function [sph, vec, coef] = SH_setup_spherical_samples(sqrt_n_samples, n_bands)
N = sqrt_n_samples;
x = rand(N,1);
y = rand(N,1);
theta = 2.0*acos(sqrt(1.0-x));
phi = 2.0*pi*y;
sph = [theta, phi, ones(length(phi),1)];
vec = [sin(theta).*cos(phi), sin(theta).*sin(phi), cos(theta)];
coef = zeros(sqrt_n_samples, n_bands^2);
for l=0:(n_bands-1)
     for m = -l:l
          index = l*(l+1)+m+1;
          coef(:,index) = SH(l,m,theta,phi);
     end
end
end

SH_project_polar_function.m

function results = SH_project_polar_function(fn, sph, coef)
weight = 4.0*pi;
fval = fn(sph(:,1), sph(:,2));
temp = repmat(fval,1, size(coef,2)).*coef;
results = sum(temp) * weight / length(sph);
end