Continuation of
http://www.math.utexas.edu/pipermail/maxima/2013/033096.html
For solving first or second order differential equations with initial
conditions
y(x0)=y0 or y(x0)=y0, y'(x0)=y1
I define function ode2_ic(eq, y, x, _ic).
Here _ic is list of initial conditions in form
_ic=[x0, y0] or _ic=[x0, y0, y1] .
(%i1)
ode2_ic(eq,y,x,_ic):=block([gsol],
gsol:ode2(eq,y,x),
if derivdegree(eq,y,x)=1 then
(
subst([x=_ic[1],y=_ic[2]],gsol),
solve(%%,%c), solve(%%,%c),
subst(%%,gsol),
logcontract(%%),
solve(%%,y)[1],
radcan(%%)
)
else
ic2(gsol, x=_ic[1], y=_ic[2],'diff(y,x)=_ic[3]))$
Example 1. y'=y*(1-y), y(0)=1/2.
(%i2) eq:'diff(y,x)=y*(1-y)$
(%i3) ode2_ic(eq,y,x,[0,1/2]);
(%o3) y=%e^x/(%e^x+1)
Example 2 . (x^2-1)*y'+2*x*y^2=0, y(0)=1.
(%i4) eq:(x^2-1)*'diff(y,x)+2*x*y^2=0$
(%i5) ode2_ic(eq,y,x,[0,1]);
(%o5) y=1/(log(x+1)+log(1-x)+1)
Example 3. y'+y/x = y^4*(1-x^2), y(1)=1.
(%i6) eq:'diff(y,x)+y/x=y^4*(1-x^2);
(%o6) 'diff(y,x,1)+y/x=(1-x^2)*y^4
(%i7) ode2_ic(eq,y,x,[1,1]);
(%o7) y=2^(1/3)/(x^(1/3)*(6*x^2*log(x)-x^2+3)^(1/3))
Example 4. 4y'' - 4y' + y = 0, y(0)=1, y'(0)=2.
(%i8) eq:4*'diff(y,x,2)-4*'diff(y,x,1)+y=0;
(%o8) 4*('diff(y,x,2))-4*('diff(y,x,1))+y=0
(%i9) ode2_ic(eq,y,x,[0,1,2]);
(%o9) y=((3*x)/2+1)*%e^(x/2)
Example 5. y'' + y = 1/cos(x), y(0)=1, y'(0)=0.
(%i10) eq:'diff(y,x,2)+y=1/cos(x);
(%o10) 'diff(y,x,2)+y=1/cos(x)
(%i11) ode2_ic(eq,y,x,[0,1,0]);
(%o11) y=cos(x)*log(cos(x))+x*sin(x)+cos(x)
Note. If solve Bernoulli equation in example 3 with ode2 and ic1 we get
wrong solution:
(%i12) eq:'diff(y,x)+y/x=y^4*(1-x^2)$
(%i13) ode2(eq,y,x)$ ic1(%, x=1, y=1);
(%o14) y=2^(1/3)/(x^(1/3)*(6*x^2*log(x)+2*x^2+3)^(1/3))
(%i15) ev(%,x=1);
(%o15) y=2^(1/3)/5^(1/3)
Correct solution is
(%i16) ode2_ic(eq,y,x,[1,1]);
(%o16) y=2^(1/3)/(x^(1/3)*(6*x^2*log(x)-x^2+3)^(1/3))
or
(%i17) expand(1/%^3);
(%o17) 1/y^3=3*x^3*log(x)-x^3/2+(3*x)/2
Aleksas D