This is straight numerical integration and provides the cdf at 1e6 equally spaced points in theta, from -pi to pi.
theta = linspace(-pi,pi,1e6);
f = exp(kappa*cos(theta-mu))/(2*pi*besseli(0,kappa));
c1 = interp1(theta,c,theta1,'spline')
c1 = 0.3148 0.7455 0.9578
For any reasonable kappa the last point of the cdf will be very close to 1, but the look of the plot depends on how far the maximum of f (at theta = mu) is from the starting point. Here the starting point is -pi, but it could be anywhere on the circle.
For an array of input angles of your choosing, interp1 provides the result.