Seemingly simple differential equation, $y'=(4x+y)/(x+4y)$ My friend has asked me for help solving the following differential equation (of which the explicit solutions are supposedly derivable):
$$\frac{dy}{dx}=\frac{4x+y}{x+4y}\tag{1}$$
I have tried hitting it with every technique I know, but none of my efforts have proved fruitful. Could you help me out? All I can end up with is a complicated implicit solution. I will now outline every approach I have taken. (Note: I have tagged every equation for referencing convenience).
[Homogeneous Substitution]
This equation is manifestly homogeneous when put in the following form:
$$\frac{dy}{dx}=\frac{4+\frac{y}{x}}{1+4\frac{y}{x}}\tag{2}$$
So let's try the substitution $u=y/x$. We can work out that $y'=xu'+u$, so plugging everything in gives us:
$$xu'+u=\frac{4+u}{1+4u}\tag{3}$$
$$\implies \frac{u'}{\frac{4+u}{1+4u}-u}=\frac{1}{x}\tag{4}$$
I can solve this by integrating both sides with respect to $x$ (left side requires heavy algebraic massaging):
$$\ln \left|(1-u)^{-5/8}(1+u)^{-3/8}\right|=\ln |x| +C\tag{5}$$
After re-substituting $u=y/x$, this then simplifies down to
$$(x-y)^5(x+y)^3=C\tag{6}$$
So the solution seems to be the solution to a eight-order polynomial, which I'm not sure can be solved for explicitly (well, it should be because why else would it show up on an entry-level DFQ homework assignment?). On top of that, a few of my previous steps have implicitly imposed domain restraints on my solution (e.g. every step where I divide by a quantity that could be zero). I've plotted the solution below, where I have taken note of the homogeneity of the original equation (the case where $C=0$ is special though - then the solutions are just $y=\pm x$).

[Integrating Factor - Exact Equation]
Ok, maybe the exact form of the solution is buried in that polynomial. Let's see if this can be made exact. Putting it in the standard form gives us:
$$f(x,y)dx+g(x,y)dy=0,$$
$$~~\textrm{where}~ f(x,y)=(4x+y)~~\textrm{and}~~ g(x,y)=-(x+4y)\tag{7}$$
This isn't an exact equation by itself ($f_y\neq g_x$). Moreover, no integrating factor of one variable (either $x$ or $y$) will work because calculating the integrating factor involves calculating an integral like:
$$\int \frac{\frac{\partial f}{\partial y}-\frac{\partial g}{\partial x}}{g(x,y)}dx\tag{8}$$
and we can't do that integral explicitly because the numerator is a constant for the specific $f(x,y)$ and $g(x,y)$ in this problem, while the denominator is a full function of $x$ and $y$. So this seems to be a no go.
[Laplace/D'Alembert Equation Form]
A [Laplace or d'Alembert equation][2] (no, not the Laplace's equation or d'Alembert's solution/formula) is first-order ordinary differential equation of the type
$$y=xf(y')+g(y')\tag{9}$$
which can be morphed into the simpler linear equation:
$$\frac{dx}{d(y')}=\left(\frac{f'(y')}{y'-f(y')}\right)x+\frac{g'(y')}{y'-f(y')}\tag{10}$$
I have found that I can manipulate the original differential equation into the desired form:
$$y=\left(\frac{y'-4}{1-4y'}\right)x\tag{11}$$
from which we can read off (well, calculate really) the corresponding differential equation for $x(y')$:
$$\frac{dx}{d(y')}=-\frac{15x}{4(1-y'^2)(1-4y')}\tag{12}$$
This is separable. Solving this gives me:
$$\ln \left|\frac{(1-y')^{5/8}(1+y')^{3/8}}{4y'-1}\right|=\ln |x| +C\tag{13}$$
this seems to be the most highly nonlinear and implicit first-order differential equation that I have ever seen. I don't know what I could realistically do from here.

~~ADDENDUM~~
These are the exact instructions from the assignment

$$\textrm{1. Find all solutions:}~~~~\frac{dy}{dx}=\frac{4x+y}{x+4y} $$

https://en.wikibooks.org/wiki/Ordinary_Differential_Equations/d%27Alembert
 A: As you note, the equation is homogeneous. The form of the coefficient matrix[*] naturally suggests introducing new variables $u$ and $v$ by
$$
\begin{aligned}
  x &= u + v, \\
  y &= u - v;
\end{aligned}\qquad
\begin{aligned}
  dx &= du + dv, \\
  dy &= du - dv.
\end{aligned}
$$
The ODE becomes
$$
\frac{du - dv}{du + dv}
  = \frac{dy}{dx}
  = \frac{4x + y}{x + 4y}
  = \frac{5u + 3v}{5u - 3v},
\tag{1}
$$
or after rearrangement,
$$
\frac{dv}{du} = -\frac{3v}{5u}.
$$
Integrating gets you quickly to $v^5u^3 = C$, the solution you found, though possibly with less effort than getting from (4) to (5) in the original post.
A level curve $f(x, y) = C$ of a smooth function satisfying an ODE is generally considered a "solution".

[*] For what it's worth, this ODE can be expressed as the linear planar system
$$
\left[\begin{array}{@{}c@{}}
    dx \\
    dy \\
  \end{array}\right]
  = \left[\begin{array}{@{}cc@{}}
    1 & 4 \\
    4 & 1 \\
  \end{array}\right]
\left[\begin{array}{@{}c@{}}
    x \\
    y \\
  \end{array}\right].
$$
The eigenvalues of the coefficient matrix are $1 \pm 4$, or $5$ and $-3$, and the vectors $(1, \pm1)$ are an eigenbasis (whose associated Cartesian coordinates are $u$ and $v$). It follows by (easy) linear algebra that parametric solutions are
\begin{align*}
\left[\begin{array}{@{}c@{}}
    x(t) \\
    y(t) \\
  \end{array}\right]
  &= \frac{1}{2}\left[\begin{array}{@{}cr@{}}
    1 &  1 \\
    1 & -1 \\
  \end{array}\right]
  \left[\begin{array}{@{}cc@{}}
      e^{5t} &  0 \\
      0 & e^{-3t} \\
  \end{array}\right]
  \left[\begin{array}{@{}cr@{}}
      1 &  1 \\
      1 & -1 \\
    \end{array}\right]
  \left[\begin{array}{@{}c@{}}
      x(0) \\
      y(0) \\
    \end{array}\right] \\
  &= \frac{1}{2} \left[\begin{array}{@{}cc@{}}
      e^{5t} + e^{-3t} & e^{5t} - e^{-3t} \\
      e^{5t} - e^{-3t} & e^{5t} + e^{-3t} \\
  \end{array}\right] \left[\begin{array}{@{}c@{}}
      x(0) \\
      y(0) \\
    \end{array}\right] \\
  &= e^{t} \left[\begin{array}{@{}cc@{}}
      \cosh(4t) & \sinh(4t) \\
      \sinh(4t) & \cosh(4t) \\
  \end{array}\right]  \left[\begin{array}{@{}c@{}}
      x(0) \\
      y(0) \\
    \end{array}\right].
\end{align*}
