Calculating lift with lifting-line theory in MATLAB
Your general formulation seems fine. Two things to point out:


If you choose to sample only half the wing, then you cannot utilize the relation: $C_L=A_1\pi A$. You would have to do the integration: $C_L=\frac{2}{V_\infty S}\int^{b/2}_{-b/2}{\Gamma(y)dy}$. That being said, maybe you can simplify the expression and get a closed-form result? Haven't tried myself.
Your code doesn't seem to handle for the case where $\sin(\theta)=0$.


You can take a look at the code below for reference:

% Wing area, span, chord etc.
S = 0.8; % [m^2], wing area
b = 1.56827; % Wing span 
s = b/2; % Semi-wing span
taper = 0.25;
Cr = S/s/(1+taper);
Ct = taper*Cr;
MAC = 2/3*Cr*(1+taper+taper^2)/(1+taper);
AR = b^2/S;

% AoA and zero-lift AoA:
AoA = deg2rad(0);
AoA_0 = deg2rad(-6.4);

% Discretization
num = 10;
[N,Ts] = meshgrid(1:num, linspace(0,pi,num));
ys = -s*cos(Ts);
Cs = (Ct-Cr)/s*abs(ys)+Cr;

% Construct linear system
B = 2*b/pi./Cs.*sin(N.*Ts)+N.*sin(N.*Ts)./sin(Ts);
if Ts(1)==0
    % Singular case, use l'Hopital's rule
    B(1,:) = N(1,:).*N(1,:);
end
RHS = AoA - ones(num,1)*AoA_0;

% Solve
A = B\RHS;
CL = A(1)*pi*AR;