dsolve('Dy = (x + y + 1)^2','x')