There is a common way of solving this kind of system translating it into a polynomial system and solving it using Gröbner basis. You can learn about it here.
First of all, as the previous answers did, I will use a different notation. I will call
$a=\cos(37º), b=\sin(37º), c=3.5\times 10^5$
$x=v_{1f}, y=v_{2f}, z=\cos\theta, t=\sin\theta$
Then, the system becomes
$
\begin{cases}
ax+yz-c=0\\
bx-yt=0\\
x^2+y^2-c^2=0
\end{cases}
$
Since $\cos\theta$ and $\sin\theta$ are related by $\cos^2\theta+\sin^2\theta=1$ we have to add the equation $z^2+t^2-1=0$, so we get
$
\begin{cases}
ax+yz-c=0\\
bx-yt=0\\
x^2+y^2-c^2=0\\
z^2+t^2-1=0
\end{cases}
$
This equations define an ideal $I=\langle ax+yz-c, bx-yt, x^2+y^2-c^2, z^2+t^2-1\rangle\subset \mathbb{R}[x,y,z,t]$, and the solutions are the elements of the algebraic set (or variety) $V(I)$.
Let's compute a Gröbner basis $G$ for $I$ using a computer algebra system (I use SAGE). Here, the monomial ordering is important, in general any lex order works fine, I will use lex with $x>y>z>t$.
$G=\{x - \frac{318367197940197676250000}{847527367068540249}t^2, y - 300000z + \frac{234667789200000}{163394518424627}t, z^2 + t^2 - 1, zt - \frac{1518598087}{2015242866}t^2, t^3 - \frac{4061203808963893956}{6367343958803953525}t\}=\{g_1,g_2,g_3,g_4,g_5\}$
Then, $V(I)=V(G)$. At first sight, the equations given by $G$ seem awful, but the last polynomial is a polynomial only on $t$, with one trivial solution $t=0$ (corresponding to $\theta=k\pi$ with $k$ integer). The other solutions are $t=\pm \sqrt{\frac{4061203808963893956}{6367343958803953525}}$ (you can clear $\theta$ using a calculator).
Then, let's go to the polynomials on $z$ and $t$. Following the notation of the book above, let's call $I_3=\langle g_5\rangle$ and $I_2=\langle g_3,g_4,g_5\rangle$. By the Extension Theorem, a partial solution $d\in V(I_3)$ extends to $(c,d)\in V(I_2)$ whenever $d\notin V(c_3,c_4)$, being $c_3,c_4$ the coefficient of the highest power of $z$ in $g_3$ and $g_4$ respectively (actually we should take $g_5$ into account, but since it has no $z$ we won't worry about it).
In this case, $V(c_3, c_4)=V(1, t)$ corresponding to $z^2$ and $z$ respectively. Since $V(1,t)=\emptyset$, all the partial solutions for $t$ that we found earlier extend to partial solutions of $V(I_2)$. So, for each solution $t$, you just have to clear $z$ in $g_3=0$ and $g_4=0$ and find common (approximate) solutions. For example, for $t=0$, $g_4=0$ for all $z$, and $g_3=0$ implies $z=\pm 1$, so here you've got two partial solutions $(1,0)$ and $(-1,0)$. Note that these values of $z$ are coherent since, coming back to your original system, $\cos\theta=\pm 1$ when $\sin\theta=0$.
Then, again, we always can extend a partial solution in $V(I_2)$ to $V(I_1)=V(\{g_2,g_3,g_4,g_5\})$ because the coefficient of $y$ (this is the maximum power of $y$) in $g_2$ is constant, so for every partial solution $(z,t)$ you can easily clear $y$ from $g_2=0$, giving a partial solution $(y,z,t)\in V(I_1)$. For example, for the partial solutions above, we've got $y=300000$ for $(1,0)$ and $y=-300000$ for $(-1,0)$, so the resulting partial solutions are $(300000, 1, 0)$ and $(-300000,-1,0)$.
Finally, extend these solutions to $V(I)$ (you can do it because the coefficient of $x$ is constant in $g_1$, and again the solutions are easy to get from $g_1=0$). Following the example above, we have found two solutions: $(0,-300000,-1,0)$ and $(0,300000,1,0)$. You can check that these are actually solutions for the polynomial system and hence from your original system. The solutions you've been given must come from one of the other possible values of $t$. In particular, from the positive value, which gives you exactly the angle $\theta=53º$.
The code to finding the Gröbner basis in SAGE, which you can try here, is the following
Q.<x,y,z,t>=PolynomialRing(QQ,order='lex')
ar = cos(37 / 360 * (2*pi)).numerical_approx()
br=sin(37 / 360 * (2*pi)).numerical_approx()
c= 3*10^5
a=ar.nearby_rational(max_error=10^(-5))
b=br.nearby_rational(max_error=10^(-5))
I=ideal(a*x+y*z-c, b*x-y*t, x^2+y^2-c^2, z^2+t^2-1 )
G=I.groebner_basis()
G
Note that these are rational approximations, which is the best we can expect from a computer, but you can make it more accurate by changing the value of max_error
.