Okay, I'll suppose this is *not* homework (i it is, please take the time to *understand* the solutions : it might be worth it...). It turns out exhibiting some interesting features (and non-features) of Sage and some packages interfaced to it, so let's dig it.

**Sage** (i. e. Maxima)

The first attempt:

```
var('x')
y = function('y')(x)
de = diff(y, x) == (4*x + y)^2
desolve(de, y, x)
```

fails with a suggestion:

```
NotImplementedError: Maxima was unable to solve this
ODE. Consider to set option contrib_ode to True.
```

which, when applied, gives a curious result:

```
Sol=desolve(de, y, ivar=x, contrib_ode=True)
Sol
[[x == _C - 1/2*arctan(1/2*sqrt(t)), y(x) == -4*x - sqrt(t)
[x == _C + 1/2*arctan(1/2*sqrt(t)), y(x) == -4*x + sqrt(t)]]
```

This parametric form introduces non only an integration constant
`_C`

but also a parameter `t`

; both `x`

and `y`

are functions of `t`

, which is only used through its
square root. We can use this to solve `sqrt(t)`

in x (hint : the justification of this solution is *not* trivial…) and get two solutions…

```
(_C, t)=var("_C, t")
SolSage=map(lambda u:u[1].subs(u[0].solve(sqrt(t))), Sol)
SolSage
```

… which turn out to be the same, and can be checked:

```
bool(SolSage[0].rhs().diff(x)==de.subs(SolSage[0]).rhs())
True
```

One notes that the constant `_C`

is determined by the initial
conditions but modulo $\pi$.

**Sympy**

By default, `Sympy`

returns a Taylor series of the solution, which
is not very interesting here. We'll ask for /all/ solutions of our
equations:

```
D=sympy.dsolve(*map(sympy.sympify, [de, y(x)]), hint="all")
D
{'1st_power_series': Eq(y(x), x**3*(C1**2*(3*C1**2 + 4) + 12*C1**2 + 16)/3 + x**5*(8*C1**2*(3*C1**2 + 8) + 4*C1**2*(9*C1**2 + 8) + C1**2*(2*C1**2*(3*C1**2 + 8) + C1**2*(9*C1**2 + 8) + 36*C1**2 + 32) + 144*C1**2 + 128)/15 + C1 + C1*x**2*(C1**2 + 4) + C1*x**4*(C1**2*(3*C1**2 + 8) + 12*C1**2 + 32)/3 + C1**2*x + O(x**6)),
'best': Eq(y(x), x**3*(C1**2*(3*C1**2 + 4) + 12*C1**2 + 16)/3 + x**5*(8*C1**2*(3*C1**2 + 8) + 4*C1**2*(9*C1**2 + 8) + C1**2*(2*C1**2*(3*C1**2 + 8) + C1**2*(9*C1**2 + 8) + 36*C1**2 + 32) + 144*C1**2 + 128)/15 + C1 + C1*x**2*(C1**2 + 4) + C1*x**4*(C1**2*(3*C1**2 + 8) + 12*C1**2 + 32)/3 + C1**2*x + O(x**6)),
'best_hint': '1st_power_series',
'default': '1st_power_series',
'lie_group': Eq(y(x), -4*x + 2*tan(C1 + 2*x)),
'order': 1}
```

The `lie_group solution can be expressed in Sage as:

```
D['lie_group']._sage_()
y(x) == -4*x + 2*tan(C1 + 2*x)
```

which is identical to the solution given by Sage up to the expression
of the constant.

**Giac**

Unfortunately, Sage is not (yet) able to translate the differential
equation in giac. But a "manual" translation allows us to call giac's
solver, which gives:

from giacpy_sage import *
SolGiac=libgiac.desolve(libgiac.diff(y(x),x)-de.rhs(),y(x))[0]
SolGiac

```
2*tan(2*x-2*c_0)-4*x
```

i. e. again the same result up to the expression of the constant…

**Fricas**

I have been unable to use Fricas solver to solve this equation. But
that's probably because I do not (yet) know the ``fricas`

' interface
in detail.

**Mathematica**

Again, Sage is currently unable to translate the differential
equation. A manual call:

```
mathematica("DSolve[Derivative[1][y][x] == (4*x + y[x])^2, y[x], x]")
{{y[x] -> -2*I - 4*x + (-I/4 + E^((4*I)*x)*C[1])^(-1)}}
```

gives a result, manually translatable to Sage as ```
-4*x +
4/(4*C1*e^(4*I*x) - I) - 2*I
```

, totally different from those obtained
so far.

A heavy pounding on this expression by deMoivre and algebraic
transformations may reveal that this expresses the same relation
between `x`

and `y`

as before. I haven't done this.

Is that homework ?

Hint :

`desolve ?`

(pay attention to the goddamn'syntax).Another hint : Sage

cangive you solutions. But not necessarily in a directly usable form, and not necessarily easy to check.