plotting autonomous differential equations

edit

Let's see (with a slight change of notations) : we start the same way you did :

sage: y=function("y")(t)
sage: DE=diff(y,t)==y^2-4*y
sage: foo=desolve(DE,y);foo
1/4*log(y(t) - 4) - 1/4*log(y(t)) == _C + t

This is only an _implicit_ solution ; not quite what we seek... Now, Sage seems reructant to transform log(a)-log(b) into log(a/b). With some reason. If we are ready to assume that, _for all t,_ y(t)!=0, we can rewrite this a bit. Let's check our substitution first :

sage: bool(foo.lhs()==1/4*(log((y(t)-4)/y(t))))
True

Okay. We can now get another equation, equivalent to the first under our assumption :

sage: bar=1/4*(log((y(t)-4)/y(t)))==foo.rhs();bar
1/4*log((y(t) - 4)/y(t)) == _C + t

Can Sage solve this ?

sage: solve(bar, y(t))
[log(y(t)) == -4*_C - 4*t + log(y(t) - 4)]

Not yet here.But the exponential has some interesting properties... Can we solve exp(bar.lhs())==exp(bar.rhs()) ?

sage: gee=bar.lhs().exp()==bar.rhs().exp();gee
((y(t) - 4)/y(t))^(1/4) == e^(_C + t)
sage: solve(gee,y(t))
[y(t) == -4/(e^(4*_C + 4*t) - 1)]

Now, we have _an_ explicit solution to your original equation. But there remain a couple of problems to solve :

• Does this describe _all_ solutions to the original equation ?

• Under which condition(s) are these solutions _real_ ?

• Are all the points of these solutions valid ?

Those questions, dear AndyH, are probably the crux of what I guess to be your _homework_, and therefore left to you as the point of the exercise...

_Computation may become easy, analysis is still hard_.

Wisdom of nations...

Here is another option for plotting. The streamline_plot does a nice job of plotting some sample curves along with arrows to give direction. I like to use the odeint solver since it works for systems so nicely. Here I use it by defining a system $ds/dt = 1, dx/dt = x^2-4x$. This gives the appropriate plot. Note that we have to tell Sage that the number 1 should be in the symbolic ring (SR) to get this to work.

var('x s t')
pp=streamline_plot(x^2-4*x,(t,0,5),(x,-4,5))
initialvalues=[-2,2]
for x0 in initialvalues:
    pp+=line(desolve_odeint([SR(1),x^2-4*x],[0,x0],srange(0,5,0.05),[s,x]),color='red',thickness=3)
pp.show(ymin=-4,ymax=5)