Is there a natural example of a discrete subgroup $\Gamma\leq PSL_2(\mathbf{R})$ such that
(1) $\Gamma$ has finite covolume
(2) $\mathfrak{h}/\Gamma$ is not compact ($\mathfrak{h}$ being the upper half-plane)
(3) $\Gamma$ is not commensurable to a conjugate of $PSL_2(\mathbf{Z})$.
I cannot think of any such example but I don't see any reason why they should not exist.
I would guess that as soon as you have more than three cusps, the Teichmuller space associated to the surface is non trivial (a complex manifold with dimension > 0) whereas there are countably many (up to isometry) surfaces whose associated group is commensurable to $PSL_2(\mathbf{Z})$. So most discrete groups associated to a surface with sufficiently many cusps (4 probably) are not commensurable to $PSL_2(\mathbf{Z})$.
All discrete subgroups of $PSL_2(\mathbf{R})$ with finite covolume are obtained as follows. Choose a 2n-gone in the hyperbolic Poincare half-plane all of whose vertices are at infinity and such that opposite sides are identified through an isometry of the half-plane. The group generated by these isometries has the polygon as its fundamental domain, and since its vertices are at infinity, this fundamental domain has finite volume. All finite covolume discrete groups without elliptic elements and with zero genus are obtained this way. This is a particular case of Poincare theorem for surfaces, that shows that there is a correspondence between finitely generated subgroups of $PSL_2(\mathbf{R})$ and convex fundamental polygons in the half-plane. See Beardon, "The Geometry of Discrete Groups", chapter 9. As soon as you have more than 6 vertices, you have a free choice for a pair of points that you cannot compensate by a global isometry (that operates triply transitively on the half-plane but no more). So you have uncountably many different subgroups that are of finite covolume and no two of them are conjugated.
EDIT: to be slightly more precise in the case of non-zero genus, all discrete groups with finite covolume are obtained by identifying sides of a finite-sided convex polygon that touches the boundary in finitely many points, and such that the angles at the vertices inside the disk that are glued together sum up to $2\pi$, or a submultiple of $2\pi$ if the vertices are fixed by elliptic elements. And all such side-paired polygons generate a discrete subgroup of finite covolume.
For explicit examples of such groups see Louis Funar's notes, page 5 (these are Hecke groups, mentioned by Anton in his comment).
Any discrete subgroup of finite covolume, whose cusps (i.e., fixed points of parabolic elements) are not the same as those of $PSL(2,{\mathbb Z})$ (i.e., ${\mathbb Q}\cup\left\{\infty\right\}$) must necessarily be not commensurable to $PSL(2,{\mathbb Z})$.
And even if the cusp set is ${\mathbb Q}\cup\left\{\infty\right\}$, the group need not always be commensurable to $PSL(2,{\mathbb Z})$. Examples of this kind have been constructed by Long and Reid in "Pseudomodular Surfaces". Their examples are non-arithmetic, which can be checked by finding elements with $tr(\gamma^2)\not\in{\mathbb Z}$.
Some more examples have been found later by Ayaka Oigo.
By playing off traces versus Dehn-Thurston coordinates, easy constructions arise.
As suggested in the comments by @YCor, a matrix representation can be written down without much trouble when expressed in terms of Dehn-Thurston coordinates on Teichmuller space that he suggested, one length-twist coordinate pair for each pants curve, except that a pants curve which is a cusp (at least one such, in order that the quotient is not compact) will not have any coordinates corresponding to it. Converting Dehn-Thurston coordinates into matrix coordinates of a corresponding representation is a straightforward one-time operation for each topological type of surface and each pants decomposition of that surface: one obtains explicit matrices for explicit generators, where the entries of the matrix are explicit functions of the lengths and twists.
As said in another comment, each group commensurate with $SL_2(\mathbb{Z})$ can be commensurated into $SL_2(\mathbb{Q})$, and therefore the traces are rational.
One also knows that for a loxodromic isometry of translation length $L$, the trace equals $2 \sinh(L/2)$ (I think I got that right...), and a closed geodesic of length $L$ corresponds to a loxodromic isometry of translation length $L$.
So just pick a number $L$ for which $2 \sinh(L/2)$ is not rational, and use $L$ as one of the pants curve lengths in the Dehn-Thurston length coordinates.
