When I run your code and store the results in a table array:
CL = 0:0.5:8;
results = table(CL.', zeros(size(CL.')), zeros(size(CL.')), 'VariableNames', {'CL', 'CR', 'solution'});
for i=1:1:17
fcn = @(CR) sind(CL(i)).^2 + sind(CR).^2 - ((0.5)*(sind(2*CL(i)))*(sind(2*CR))) - (2*(sind(CL(i)).^2)*(sind(CR).^2)) - sind(CL(i)-CR).^2;
CR(i) = fzero(fcn,1 );
%disp(CR);
results{i, 'CL'} = CL(i);
results{i, 'CR'} = CR(i);
results{i, 'solution'} = fcn(CL(i));
end
displaying that table shows that the value of your function for the solution found by fzero is extremely small for each value of CL.
results =
17×3 table
CL CR solution
___ _______ ___________
0 1 0
0.5 1 2.3366e-20
1 0.97978 -4.124e-20
1.5 0.97376 -2.6978e-19
2 0.97172 -3.3119e-19
2.5 1.04 7.5217e-19
3 0.97172 8.0976e-19
3.5 0.9755 8.3687e-19
4 1 1.0639e-18
4.5 0.98637 -6.0986e-19
5 0.97172 -2.575e-19
5.5 -1.2817 -2.3852e-18
6 -2.5999 -2.4937e-18
6.5 0.96 -2.6563e-18
7 0.99967 -6.2884e-18
7.5 -6.0998 -2.9273e-18
8 0.98114 -7.5894e-19
The values of CR found by fzero may not be the solution you expected, but they seem from these results to each be a solution to your equation for the given value of CL.
If we plot one of your functions, I see one reason why you may be confused. I'm using just one of the CL values over which you iterated, so I was able to simplify your function by removing the indexing on CL:
CL = 4;
fcn = @(CR) sind(CL).^2 + sind(CR).^2 - ((0.5)*(sind(2*CL))*(sind(2*CR))) - (2*(sind(CL).^2)*(sind(CR).^2)) - sind(CL-CR).^2;
CR2 = linspace(-3, 3, 1000);
plot(CR2, fcn(CR2))
Look at the Y limits. The maximum value of your function over that interval is about 7e-18. That's pretty small.
