A particularly simple non-homogeneous example in which one can explicitly integrate the Jacobi equations is the complete metric on $\mathbb{R}^2$ given by $$ g = (x^2{+}y^2{+}2)\bigl(\mathrm{d}x^2+\mathrm{d}y^2\bigr). $$ It has Gauss curvature $K = -4/(x^2{+}y^2{+}2)^3<0$, and, visibly, a rotational symmetry about the origin $(x,y)=(0,0)$.
It is not hard to show that, up to a rotation, each geodesic can be parametrized in the form $$ (x,y) = \bigl(r\,\cosh t,\ \sqrt{r^2+2}\,\sinh t\,\bigr) $$ where the constant $r\ge0$ determines the closest approach of the geodesic to the origin. The element of arc length along this geodesic is then found to be $\mathrm{d}s$, where $$ s = t + (r^2{+}1)\,\cosh t\,\sinh t. $$
Now, the Jacobi fields split into the tangential Jacobi fields, which are spanned by $$ J_1 = \frac{\partial}{\partial s} =\frac{1}{(1+(r^2{+}1)\cosh 2t)}\,\frac{\partial}{\partial t} \quad\text{and}\quad J_2 = s\,\frac{\partial}{\partial s}, $$ and the normal Jacobi fields $J_3= f_1\,N$ and $J_4 = f_2\,N$, where $N$ is the unit normal vector field to the curve and $f_1$ and $f_2$ are a basis for the solutions to the (linear) normal Jacobi equation $$ \frac{d}{ds}\left(\frac{df}{ds}\right) + K\,f = 0. $$ Using the above formuale, one finds that these can be taken to be $$ f_1(t) = r^2+1+\cosh 2t\quad\text{and}\quad f_2(t) = \sinh 2t\,. $$
Finally, note that these formulae generalize immediately to the case of the cohomogeneity-1 metric on $\mathbb{R}^n$ with the formula $$ g = \bigl(|x|^2+2\bigr)\,(\mathrm{d}x\cdot\mathrm{d}x), $$ since every geodesic in this space lies in a 2-plane through the origin $x=0$. In fact, the more general family of complete, conformally flat Liouville metrics $$ g = \bigl(a_0 + a_1\,{x_1}^2 + \cdots + a_n\,{x_n}^2\bigr)(\mathrm{d}{x_1}^2 + \cdots + \mathrm{d}{x_n}^2), $$ where $a_0>0$ and $a_i\ge 0)$ for $1\le i\le n$, has the property that the Jacobi equations on the geodesics of this metric can be explicitly integrated.
Remark: For more information about the metric $g$, in particular the explicit formula for its distance function, see this answer of mine
Additional Examples: In case the OP is interested, here is another group of examples that may be of interest. These are also conformally flat Liouville metrics, but now defined on the interior of the unit $n$-ball $B = \{\,x\in\mathbb{R}^n\ |\ |x|\le 1\ \}$ for constants $m_i>0$ $(1\le i\le n)$ as metrics of the form $$ g = \left(1-|x|^2\right)\left(\frac{{\mathrm{d}x_1}^2}{{m_1}^2}+\cdots + \frac{{\mathrm{d}x_n}^2}{{m_n}^2}\right) $$ It is not hard to show that every $g$-geodesic is parametrized in the form $$ x_i(t) = \lambda_i\,\cos( m_i\,t + q_i) $$ where the constants $\lambda_i$ and $q_i$ satisfy ${\lambda_1}^2+\cdots+{\lambda_n}^2 = 1$ and $q_1+\cdots+q_n=0$ and where arclength $s$ along the geodesic satisfies $$ \mathrm{d}(vs+c) = \bigl(1-|x(t)|^2\bigr)\,\mathrm{d}t $$ for some constants $v$ and $c$. Note that the constants $(\lambda,q,v,c)$ subject to the two constraints vary in a manifold of dimension $2n$, which is the dimension of the space of (parametrized) $g$-geodesics in $B$. Now, the Jacobi fields along a given geodesic can be computed explicitly as the partials of the explicit formulae with respect to the constrained parameters $(\lambda,q,v,c)$ using the Chain Rule and eliminating $s$ in favor of $t$, which can be done explicitly using the given relation between $\mathrm{d}s$ and $\mathrm{d}t$.