A doubt on calculation of tetrad basis vectos of a non-diagonal metric tensor First of all, I will ask the patience of community because this is a "need-an-explicit-calculation answer". Pure abstract considerations will not, I think, help me very much.
So, I would like some sort of step-by-step rule to calculate tetrad coefficients. I will explain (please consider the metric tensor signature as $(-,+,+,+)$:
PART I
For diagonal metric tensors, we have that the tetrad basis vectos are given by:

$$\begin{cases}\hat{\textbf{e}}_{(0)} = \frac{1}{\sqrt{-g_{00}}}\frac{\partial}{\partial x^{0}}\\ \hat{\textbf{e}}_{(1)} = \frac{1}{\sqrt{g_{11}}}\frac{\partial}{\partial x^{1}}\\\hat{\textbf{e}}_{(2)} = \frac{1}{\sqrt{g_{22}}}\frac{\partial}{\partial x^{2}}\\\hat{\textbf{e}}_{(3)} = \frac{1}{\sqrt{g_{33}}}\frac{\partial}{\partial x^{3}}\end{cases} \tag{1}$$

So concerning an "engieneer way of think", we can say that:

Given a (diagonal) metric tensor, apply the formulas $(1)$, and then you receive as output the tetrad basis.

PART II
So, I would like to build  a some sort of a general algorithm (a.k.a. some calculation which gives us the right results, always), as in PART I, but now for a non-diagonal metric tensor. The metric tensor is given by $[1]$:

$$ds^2 = -N(r,\theta)^2dt^2 + \frac{1}{1-\frac{b(r,\theta)}{r}}dr^2+r^2K^2(r,\theta)\Bigg\{d\theta^2+sin^2(\theta)\bigg[d\phi-\omega(r,\theta)dt\bigg]^2\Bigg\} $$

And the tetrad basis of the paper $[2]$ is given by:

\begin{cases}\hat{\textbf{e}}_{(0)} = \frac{1}{N}\frac{\partial}{\partial x^{0}} + \frac{\omega}{N}\frac{\partial}{\partial x^{3}}\\ \hat{\textbf{e}}_{(1)} = \sqrt{1-\frac{b(r,\theta)}{r}}\frac{\partial}{\partial x^{1}} \\ \hat{\textbf{e}}_{(2)} = \frac{1}{rK}\frac{\partial}{\partial x^{2}} \\ \hat{\textbf{e}}_{(3)} = \frac{1}{rKsin(\theta)}\frac{\partial}{\partial x^{3}} \tag{2} \end{cases}

PART III
So my doubt begins when I say that, I have no idea on how to obtain the vectors in $(2)$. Also I do not know if it is possible to stablish an general way to calculate tetrad vectors given any type of metric tensor. I have applied, in Mathematica, the function
$$MatrixForm[FullSimplify[DiagonalMatrix[Eigenvalues[metric]]]]$$
But I didn't be able to go any further. So, my doubt is:

How can I calculate the vectors in (2)?

I appreciate an step-by-step answer, as I said above, but please is not mandatory.
$$ * * * $$
$[1]$ TEO.E. Rotating Traversable Wormholes https://arxiv.org/abs/gr-qc/9803098
$[2]$ LOBO.F.S.N. Exotic Solutions in General Relativity https://arxiv.org/abs/0710.4474
 A: I believe that for ANY given metric tensor (diagonal or not) the vector you give in (1) will be the elements of the base (they are also normalized) in which your metric tensor is written but they are just Linearly Independent, but no necessary orthogonal.
They will only be orthogonal if the metric tensor you are given is diagonal, otherwise, like in the Kerr metric you've written you can see that in (1) e_0 is not orthogonal to e_3 utilizing the Kerr metric in the inner product.
Now, the vectors in (2) are just the gram schmidt algorithm to create an orthonormal basis (2) from the non orthogonal basis (1).
And now if you write the kerr metric on the new basis (2) you will see that the matrix will be diagonal.
An easier way to see this is to make the substitution $t=Nt'-w\phi$ and see that there will be cross terms between $dt'$ and $d\phi$
A: The standard way to get an orthonormal basis (tetrad) from an arbitrary linearly independent basis (e.g. coordinate basis) in an inner product (metric) space is Gram-Schmidt orthonormalization.
You can simply do this procedure at each point, starting from the coordinate basis, giving you one of the many possible orthonormal tetrads. If you want a different one you can apply a local Lorentz transformation.
A: This is how I would do it. Let my $g$ be my metric tensor. With respect to a cooridnate system of choice, I can write this down as
$$ g = g_{\mu \nu} \mathrm{d}x^\mu \otimes \mathrm{d}x^\nu$$
as you have in part II of your question, where $\mu , \nu$ are coordinate indices. Define the tetrad basis vectors as
$$ e_a = e_a^{\ \mu} \frac{\partial}{\partial x^\mu}, \quad a = 0,1,2,3. $$
By definition, the tetrad is an orthonormal basis so they must satisfy $g(e_a , e_b) = \eta_{ab}$, where $\eta_{ab} = \mathrm{diag}(-1,1,1,1)$ is the Minkowksi metric. In terms of components, this relations reads
$$ e_a^{\ \mu} e_b^{\ \nu} g_{\mu \nu} = \eta_{ab}$$
In other words, the tetrad basis diagonalises the metric. So finding the tetrad basis is equivalent to finding a matrix $e_a^{\ \mu}$ which diagonalises the matrix $g_{\mu \nu}$! Once you have found this, you can plug the components $e_a^{\ \mu}$ back into the expression for the tetrad above.
For your example, you have a metric of the form
$$ g = -A^2 \mathrm{d}t^2 + B^2 \mathrm{d}r^2 + C^2 \left[ \mathrm{d}\theta^2 + D^2(\mathrm{d}\phi - E \mathrm{d}t )^2 \right] $$
The components $g_{\mu \nu}$ in your case are
$$ g_{\mu \nu} = \begin{pmatrix} -A^2 + C^2 D^2 E^2 & 0 & 0 & -C^2D^2E \\ 0 & B^2 & 0 & 0 \\ 0 & 0 & C^2 & 0 \\ -C^2 D^2 E & 0 & 0 &C^2 D^2  \end{pmatrix}   $$
which is diagonalised by
$$ e_a^{\ \mu} = \begin{pmatrix} 
\frac{1}{A} & 0 & 0 & \frac{E}{A} \\
 0 & \frac{1}{B} & 0 & 0 \\ 
0 & 0 & \frac{1}{C} & 0 \\
 0 & 0 & 0 & \frac{1}{CD} 
\end{pmatrix} $$
which gives you the components of your four tetrad basis vectors in your equation (2).
As metric tensors are symmetric tensors, you can always diagonalise them to find a tetrad and the components of the tetrad are given by the rows of the matrix that diagonalises it.
By the way, in general you won't be able to find a single tetrad field that covers your entire spacetime, otherwise that would signify your spacetime is actually trivally Minkowski! Also, the tetrad basis not unique and is defined up to a Lorentz transformation, i.e. if $\{ e_a \}$ is a tetrad basis, then so is $\{ e'_a = \Lambda_a^{\ b} e_b \}$, where $\Lambda $ is a Lorentz transformation.
