You are essentially there. Wien's law constricts the form of the blackbody spectrum to $$\rho(\lambda,T)=\frac{f(\lambda T)}{\lambda^5},$$ while Wien's displacement law talks about the peak of the spectrum at a given temperature. Thus you are correct that you need to set $\cfrac{\partial \rho}{\partial\lambda}=0$, and that this leads to the equation
$$5 f(\lambda T)=\lambda Tf'(\lambda T).\tag 1$$
This equation will look a lot more friendly if you introduce some more notation. If you set $\mu=\lambda T$, then you can reduce (1) to an equation that's exclusively in $\mu$:
$$5 f(\mu)= \mu f'(\mu).\tag 2$$
Depending on $f$, this may have one, many, or no solutions, and each solution will mark a peak in the blackbody spectrum. Since Wien's law does not specify $f$, we can't tell yet, but we expect that there will be a unique solution $\mu_0$ to (2). This implies that the peak wavelength $\lambda_\text{peak}$ must obey
$$\lambda_\text{peak}=\frac{\mu_0}T.$$
Note that here $\mu_0$ is a (dimensionful) constant that is determined by the final form of $f$, but we know it can't depend on anything other than fundamental constants.
It's also important to note that this fact is part of a stronger result which is what Wien's law really embodies: a change in temperature can only make a scaling transformation on the corresponding blackbody spectrum. This is beautifully explored in Using Wien's Law to show spectral distruibution function of one temperature represents all temperatures.
But this is for perfect blackbodies only, which have no theoretical existence.
Just for clarification, the Wien displacement law is not valid for emission spectrum "of blackbodies only". It is valid for any radiation in thermodynamic equilibrium with matter (equilibrium radiation), which is a more general concept that does not require presence of any black body. This equilibrium radiation has energy density spectrum given by the Planck function
$$
\rho_{Planck}(\lambda, T) = \frac{8\pi h c}{\lambda^5}\frac{1}{e^{\frac{hc}{\lambda kT}}-1}
$$
and can be produced and maintained by non-black bodies, for example by enclosure of any material brought to thermodynamic equilibrium at temperature $T$. Wien's displacement law says that the product of $T$ and $\lambda_{max}$ for which the energy density spectrum is maximum is the same for all temperatures.
Does a similar formula exist for real bodies, which expresses λT in terms of its emissitivity ϵ?
Let me rephrase your question : for emission intensity spectrum $I(\lambda, T)$ of non-black body with known emissivity $\epsilon(\lambda, T)$, is there a simple relation between the wavelength of maximum and temperature, similar to above Wien's displacement law ?
The answer is that in general, most probably no, because the emission radiation of real body has intensity spectrum proportional to
$$
I(\lambda, T) \propto \rho_{Planck}(\lambda, T) \epsilon(\lambda, T),
$$
which means that the number, positions and height of peaks are modified by the emissivity $\epsilon(\lambda, T)$. This depends on the body, its material and thickness and can be fairly general function of $\lambda$ and $T$. There may be some general restriction on this function due to physical laws, but Planck's formula does not imply any, apart from the condition that $0<\epsilon \leq 1$ for all frequencies. One would have to study more detailed model of matter explaining value of $\epsilon(\lambda, T)$ to find such relations.
Best Answer
Link to original version (German)
Willy Wien: Über die Energievertheilung im Emissionsspectrum eines schwarzen Körpers. Annalen der Physik, Band 294, Nr. 8, S. 662–669 (1896)
http://myweb.rz.uni-augsburg.de/~eckern/adp/history/historic-papers/1896_294_662-669.pdf
Today's derivation
We calculate the eigenmodes of a box, where the mode index is $j^2 = j_x^2+j_y^2+j_z^2=\left(\frac{2\nu}{c}L\right)^2$ where we used the condition of resonance.
We calculate the number of modes $G(\nu)=2\frac{1}{8}\frac{4\pi}{3}j^3$ in the frequency spectrum between 0 and \nu.
We calculate the spectral mode density $g(\nu)=\frac{\partial G(\nu)}{\partial\nu}$. The spectral energy density $u(\nu)$ is now the product of $g(\nu)$ and the energy per mode $\epsilon_{Wien}=h\nu e^{-\frac{h\nu}{k_B T}}$ (from classical Boltzmann statistics) per volume $L^3$.
The rest (simple math) is up to you or check the German reference 1. You get the relation from above and $$u(\nu)=\rho(\nu)e^{-\frac{h\nu}{k_B T}}\;.$$