> # This is an example of a DE which Maple will solve but the solution is composed of # Bessel functions. However, plot can be used to get an idea of what the solution # looks like. Also, DEplot is used as a comparison. Finally, numeric dsolve is used to # show an answer accurate to 7 decimal places. # # Written by Doug Hensley # # Solve y' = x^2 - y^2 # > dsolve({diff(y(x),x)=x^2-y(x)^2,y(0)=1},y(x)); 2 2 1/2 2 / y(x) = 2 x GAMMA(3/4) (GAMMA(3/4) 2 - Pi) BesselK(-3/4, 1/2 x ) / (( / 2 2 1/2 2 GAMMA(3/4) (GAMMA(3/4) 2 - Pi) BesselK(1/4, 1/2 x ) - 2 -------------------------------------------------------- 4 2 Pi (- 2 GAMMA(3/4) + Pi ) 2 4 2 + BesselI(1/4, 1/2 x )) Pi (- 2 GAMMA(3/4) + Pi )) + 2 / x BesselI(-3/4, 1/2 x ) / ( / 2 2 1/2 2 GAMMA(3/4) (GAMMA(3/4) 2 - Pi) BesselK(1/4, 1/2 x ) - 2 -------------------------------------------------------- 4 2 Pi (- 2 GAMMA(3/4) + Pi ) 2 + BesselI(1/4, 1/2 x )) -------------------------------------------------------------------------------- > expn:=simplify(rhs("));\ 4 2 1/2 expn := (2 GAMMA(3/4) BesselK(-3/4, 1/2 x ) 2 2 2 - 2 GAMMA(3/4) BesselK(-3/4, 1/2 x ) Pi 2 4 3 2 - 2 Pi BesselI(-3/4, 1/2 x ) GAMMA(3/4) + Pi BesselI(-3/4, 1/2 x )) x / 4 2 1/2 / (- 2 GAMMA(3/4) BesselK(1/4, 1/2 x ) 2 / 2 2 + 2 GAMMA(3/4) BesselK(1/4, 1/2 x ) Pi 2 4 2 3 - 2 BesselI(1/4, 1/2 x ) Pi GAMMA(3/4) + BesselI(1/4, 1/2 x ) Pi ) -------------------------------------------------------------------------------- > ?BesselK -------------------------------------------------------------------------------- > plot(expn,x=-2..2); > -------------------------------------------------------------------------------- > expn2:=map(evalf,expn); 2 expn2 := ( - 3.057184416 BesselK(-.7500000000, .5000000000 x ) 2 / + 16.83806459 BesselI(-.7500000000, .5000000000 x )) x / ( / 2 3.057184416 BesselK(.2500000000, .5000000000 x ) 2 + 16.83806459 BesselI(.2500000000, .5000000000 x )) -------------------------------------------------------------------------------- > evalf(subs(x=2,expn2)); 1.679458983 -------------------------------------------------------------------------------- > -------------------------------------------------------------------------------- > with(DEtools); [DEplot, DEplot1, DEplot2, Dchangevar, PDEplot, dfieldplot, phaseportrait] -------------------------------------------------------------------------------- > DEplot1(x^2-y^2,[x,y],-2..2,{[0,-2],[0,-1],[0,-1/2],[0,0],[0,1],[0,2]},-2..2,limitrange=true); > -------------------------------------------------------------------------------- > ?DEtools[options] -------------------------------------------------------------------------------- > Digits:=12;evalf(subs(x=2,expn2)); -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- > fn:=dsolve({diff(y(x),x)=x^2-y(x)^2,y(0)=1},y(x),numeric); fn := proc(rkf45_x) ... end -------------------------------------------------------------------------------- > fn(2); [x = 2, y(x) = 1.679458972594270] -------------------------------------------------------------------------------- >