Hello,

Consider the following differential equation

```
y=function('y')(x)
eq=6*diff(y,x)-2*y==x*y^4
sol=desolve(eq,y,show_method=True)
```

This code gives the following solution

```
[e^(1/3*x)/(-1/2*(x - 1)*e^x + _C)^(1/3), 'bernoulli']
```

When I solve this with Mathematica it gives me 3 solutions for y(x) (for each root, as expected).

Question #1: Why Sage gives one solution? When I solve this differential equation analytically, I would multiply both sides by y^(-4) (assuming y!=0) and then change the variable y^3=v. The differential equation becomes a 1st order linear differential equation in the dependent variable v and independent variable x. After solving for v, the cube-root brings the 3 solutions Mathematica is giving.

Question #2: Ok, this is really bugging me the most: Consider the same differential equation but this time let's also give the initial condition y(0)=-2

```
y=function('y')(x)
eq=6*diff(y,x)-2*y==x*y^4
sol=desolve(eq,y,show_method=True,ics=[0,-2])
```

This gives the following

```
[e^(1/3*x)/(-1/2*(x - 1)*e^x - 1/2)^(1/3), 'bernoulli']
```

However, this function blows up at x=0! How can I have a solution by calling desolve with the initial condition y(0)=-2 and end up with a solution that doesn't satisfy the initial condition ?!

Question #3: I tried to figure out the constant by myself with the following code:

```
y=function('y')(x)
eq=6*diff(y,x)-2*y==x*y^4
sol=desolve(eq,y,show_method=True)
c=sol[0].variables()[0] #<--- this refs _C
icsol=solve(sol[0].subs(x==0)==-2,c)
```

icsol gives me

```
[(_C + 1/2)^(1/3) == (-1/2)]
```

Why doesn't it give as [_C==...] ?

Question #4: For the code given in Question #3, if I try to substitute icsol back in the solution as

```
sol[0].subs(icsol)
```

I got sol[0] thrown back at me unchanged since there is no (_C + 1/2)^(1/3) in sol[0].

How can I even substitute back the constant C in the solution?

Thanks in advance for your help.