Why I am trying what I am trying in my ode program

> First there is probably a way to symbolically derive equations for the terms of Taylor series terms from the differential equations in terms of the initial conditions. 

For one first order DE, try something like:

(%i1) taylor_solve_de(de, x, y, [l]) := block([xo,yo,n,k],
  de : lhs(de)-rhs(de),
   xo : assoc('xo,l),
   yo : assoc('yo,l),
   n : assoc('n,l),
   taylor_solve_de_xx(de,y,x,yo,xo,1,n))$

(%i2) taylor_solve_de_xx(de,y,x,yo,xo,k,n) := block([g,eq,sol],
   if k >= n then yo else (
       g : gensym("g"),
       eq : subst(y = yo + g * (x-xo)^k,de),
       eq : ev(eq,diff),       
       eq : tlimit(eq / (x-xo)^(k-1),x,xo),
       sol : algsys([eq],[g]),
       taylor_solve_de_xx(de, y, x, yo + subst(sol, (x-xo)^k*g),xo,k+1,n)))$

Taylor solution to y' = mu * y:

  (%i3) taylor_solve_de('diff(y,x) - mu * y, x, y, xo=0,yo=1,n=4);
  (%o3) (mu^3*x^3)/6+(mu^2*x^2)/2+mu*x+1

Order nine solution to ' diff(y,x) = x*(1+y)^(1/2), x,  with initial condition (x=0,y=12)

  (%i4) taylor_solve_de('diff(y,x) = x*(1+y)^(1/2), x, y, xo=0,yo=12,n=9);
  (%o4) x^4/16+(sqrt(13)*x^2)/2+12

What happened to those O(x^5) terms?  Let's try order 19:

  (%i5) taylor_solve_de('diff(y,x) = x*(1+y)^(1/2), x, y, xo=0,yo=12,n=19);
  (%o5) x^4/16+(sqrt(13)*x^2)/2+12

Same approximation....oh, the DE has a solution that is a 4th degree polynomial :)

  (%i6) subst(y = %, 'diff(y,x) - x*(1+y)^(1/2));
  (%o6) 'diff((x^4/16+(sqrt(13)*x^2)/2+12),x,1)-x*sqrt(x^4/16+(sqrt(13)*x^2)/2+13)

  (%i7) %,diff;
  (%o7) -x*sqrt(x^4/16+(sqrt(13)*x^2)/2+13)+x^3/4+sqrt(13)*x

  (%i8) scanmap('gfactor,%);
  (%o8) 0

(1) I thought I could generate the Horner form automatically, but I didn't get that to work.  I also think the code 
     could be modified to generate continued fraction solutions.  Variations are possible too (Frobenius-like solutions).

(2)  Of course, this could be boosted up to systems of first order DEs.

(3) This code will work poorly with floating point numbers---maybe the intent of your code was to do things like
     this, but stick purely to floating point calculations (for speed)?

(4) Likely my code is missing a few things that I can't think of right now.

--Barton