N = 4;

%set up matrix:
M = zeros(length(x),1+2*N);
M(:,1) = 1;
for k=1:N
    M(:,2*k) = cos(2*pi/L*k*x);
    M(:,2*k+1) = sin(2*pi/L*k*x);
end

C = M\y; %equivalent to (M'*M)\(M'*y)

% least squares coefficients:
A0L2 = C(1);
AL2 = C(2:2:end);
BL2 = C(3:2:end);

%reconstruct the function values
yyL2=A0L2;
for k=1:length(AL2)
    yyL2 = yyL2 + AL2(k)*cos(2*pi/L*k*x) + BL2(k)*sin(2*pi/L*k*x);
end