just implement the Modified Bessel function. it's easy.
i always like my window definitions centered about zero, since pretty much all of them are even symmetry.
i'll do this in discrete-time, but it's essentially the same thing in continuous-time.
Kaiser window:
$$ w[n] \triangleq \begin{cases} \frac{1}{I_0(\beta)} I_0\left(\beta \sqrt{1 - \left(\frac{n}{M/2}\right)^2} \right) \quad \quad & |n| \le M/2 \\
0 & |n|>M/2 \\
\end{cases} $$
$I_0(x)$ is the 0th-order modified Bessel function of the 1st kind. $M+1$ is the number of non-zero samples or FIR taps (the FIR filter order is $M$ and, in my centered and symmetrical case, must be even). $\beta$ is a "shape parameter" and O&S recommend this heuristic:
$$ \beta = \begin{cases} 0.1102 \cdot (A-8.7) & A>50 \\
0.5842 \cdot (A-21)^{2/5} + 0.07886 \cdot (A-21) \quad & 21 \le A \le 50 \\
0.0 & A<21 \\
\end{cases}$$
$$ M = 2 \left\lceil \frac{A-8}{4.57 \cdot \Delta\omega} \right\rceil $$
$A$ is the desired stopband attention in dB and $\Delta\omega$ is the desired width of the transition band in normalized angular frequency.
finally, the Bessel is evaluated as:
$$ I_0(x) = 1 \ + \ \lim_{K \to \infty} \ \sum\limits_{k=1}^{K} \left(\frac{x^2}{4}\right)^{k} (k!)^{-2} $$
when you evaluate this with a computer, pick a $K$ decently large (my guess is that $K=32$ is good enough) and evaluate the summation starting with $k=K$ and work it backwards to $k=1$ to keep numerical accuracy. you might want to use Horner's method.
$$ I_0(x) \approx 1 + x^2\left( \tfrac{1}{(1!)^2 \, 4^1} + x^2\left(\tfrac{1}{(2!)^2 \, 4^2} + x^2\left(... + \, x^2\left(\tfrac{1}{((K-1)!)^2 \, 4^{K-1}} + x^2 \tfrac{1}{(K!)^2 \, 4^K} \right) \right) \right) \right) $$
you can evaluate all of the $(k!)^{-2}$ in advance with a short table.
Alright, someone made me do some work. This took about 45 minutes to code up and debug. So here is my MATLAB code for implementing the 0th-order Bessel function of the first find (which is $I_0(x)$ above):
function y = mybessel(x)
%
% Computes the 0th-order Modified Bessel function of the first kind
%
K = 32;
bessel_coef = zeros(1,K);
kfac = 1;
two_to_the_k = 1;
for k = 1:K
kfac = kfac * k;
two_to_the_k = two_to_the_k * 2;
bessel_coef(k) = 1/(kfac*two_to_the_k)^2; % compute power series coefficients in advance
end
x = x.^2;
y = x .* bessel_coef(K);
for k = K-1:-1:1
y = x .* (bessel_coef(k) + y); % Horner's method
end
y = 1 + y;
end
and here is the test code:
x = linspace(-16, 16, 32*4096+1);
I_0 = besseli(0, x); % MATLAB's modified bessel
y = mybessel(x); % my bessel
figure(1)
plot(x, I_0)
hold on
plot(x, y)
hold off
figure(2)
plot(x, y - I_0) % plot error
with results:

and error:

for $|x| \le 16$, the relative error is less than $10^{-15}$. the error increases with increasing $|x|$.