Solution 1:

Separation of variables: let $f(u,v)=U(u)V(v)$. Substituting into the PDE and simplifying we get

$$ \frac{1}{U}\left[uU'+u^2 U''-U'' \right]=-\frac{1}{V}\left[vV'+v^2 V''-V'' \right] $$

Introducing the arbitrary constant $k^2$

$$\tag{1} uU'+(u^2-1)U''=k^2U $$ $$\tag{2} vV'+(v^2-1)V''=-k^2V $$

Let us study (1). The results will carry over to (2) with the replacement $k^2 \to-k^2$. This ODE may be cast in Sturm–Liouville form

$$\tag{3} \frac{d}{du}\left( \sqrt{1-u^2} U'\right) =-k^2 \frac{U}{\sqrt{1-u^2}} $$

With weight $w(u)=1/\sqrt{1-u^2}$ and eigenvalue $k^2$. For a boundary value problem posed on $-1<u<1$, eq (3) is a regular Sturm–Liouville equation$^\dagger$. With the assistance of Mathematica, we find the two solutions

$$\tag{4} U^1_k(u)=\sin \left[k \tan^{-1} \left( \frac{u}{\sqrt{1-u^2}}\right)\right] \\ U^2_k(u)=\cos\left[k \tan^{-1} \left( \frac{u}{\sqrt{1-u^2}}\right)\right] $$

The general solution to (1) is then

$$\tag{5} U(u)=\sum\limits_k (A^1_k U^1_k(u) + A^2_k U^2_k(u)) $$

With coefficients $A^1$ and $A^2$ chosen to match boundary conditions. The sum is over all allowed values of $k$ consistent with the boundary data. For example, with the Dirichlet problem $U(0)=0$, $U(1)=0$, all $A^2$ vanish, and $k=2n$, for $n \in \mathbb{Z}$

We can find the solutions to (2) with the replacement $k\to ik$ in (4)

$$\tag{6} V^1_k(v)=\sinh \left[k \tan^{-1} \left( \frac{v}{\sqrt{1-v^2}}\right)\right] \\ V^2_k(v)=\cosh\left[k \tan^{-1} \left( \frac{v}{\sqrt{1-v^2}}\right)\right] $$

The general solution for $f$ is

$$\tag{7} f(u,v)=\sum\limits_k (A^1_k U^1_k(u) + A^2_k U^2_k(u))(B^1_k V^1_k(v) + B^2_k V^2_k(v)) $$

In order to actually find coefficients once boundary data is set, you will need to use the orthogonality of the eigenfunctions with weight $w$. For the $U$ equation, these read

$$\tag{8} \int\limits_{u_1}^{u_2}du\ (1-u^2)^{-1/2} U^i_k(u)U^j_{k'}(u) =\delta^{ij}\delta_{k,k'} $$

Where $i,j=1,2$, and $u_1<u<u_2$. An identical equation follows for the $V$. Continuing the example, if $V(0)=0$ and $V(1)=1$ we may apply (8) to find $A^1_n = \frac{\sin^2(n \pi /2)}{n \sinh(n \pi)}$. Thus the solution with these boundary conditions and $U(0)=0$, $U(1)=0$ is

$$\tag{9} f(u,v)=\sum\limits_{n} \frac{ \sin^2(n \pi/2)}{n \sinh(n\pi)}\sin(2n \tan^{-1}(u/\sqrt{1-u^2}))\sinh(2n \tan^{-1}(v/\sqrt{1-v^2})) $$

Here is a plot of the solution in (9)

enter image description here

$\dagger$ For $u>1$ the equivalent of (3) reads

$$ \frac{d}{du}\left( \sqrt{u^2-1} U'\right) =k^2 \frac{U}{\sqrt{u^2-1}} $$

In this case, the eigenvalue is $-k^2$, the weight $w(u)=1/\sqrt{u^2-1}$, and in the solutions (4), all trig/ hyperbolic functions are to be exchanged.