Numerically you won't be able to find the exact values, more so as the data has 2 digit precision and the values you want to obtain have 4 digit precision.
As you can see below, for tolerance of 0.1 for step-size 0.01, you get 2 out of 3 intersection points -
a=0:0.01:9;
y1=sqrt(64-a.^2);
y2=-1.*a.*cot(a);
a(abs(y1-y2)<0.1)
For step-size 1e-4 and tolerance 1e-3, you get 4 different roots
a=0:0.0001:9;
y1=sqrt(64-a.^2);
y2=-1.*a.*cot(a);
a(abs(y1-y2)<0.001)
Even solve() is unable to find a solution symbolically, and returns a numeric solution via vpasolve. Note that vpasolve only returns a single solution for non-polynomial equations
syms x
assume(x>=0 & x<=9)
f1(x)=sqrt(64-x.^2);
f2(x)=-1.*x.*cot(x);
sol = solve(f1-f2==0, x)
You can call vpasolve() multiple times by setting 'Random' as true to obtain the solutions -
for k=1:5
s = double(vpasolve(f1-f2==0, x, [0 9], 'Random', true))
end
plot(a,y1)
grid on
hold on
plot(a, y2)
ylim([0, 12])
