Code similar to the commercial macsyma's CONVERT_ODES_TO_FIRST_ORDER_SYSTEM



Hello All,
Following is a function to convert ode's to a first order system,
and an example from the commercial Macsyma's documentation.
As usual, I am interested in any comments.
Dan Stanger
convert_to_first_order_system(exp,deplist,indep,[extra]):=block(
   [subs:[],diffs:[]],
   if not listp(exp) then exp:[exp],
   if not listp(deplist) then deplist:[deplist],
   if length(extra) = 1 then ( /* keep consistant with the commercial 
product */
      if listp(first(extra)) then extra:first(extra)),
   for d in deplist do block([dd],
     push('diff(d,indep),diffs),
     dd:derivdegree(exp,d,indep),
     for di:2 thru dd do block([nv,s1,s2],
        nv:concat(d,"_",di-1),
        subs:push(nv='diff(d,indep),subs),
        s1:'diff(d,indep)=nv,

        exp:subst(s1,exp),
        extra:subst(s1,extra),
        s2:'diff(d,indep,di)='diff(nv,indep,di-1),
        push(rhs(s2),diffs),
        exp:subst(s2,exp)
        )),
   [solve(append(exp,subs),diffs),subs,extra]
);

des:['diff(chi,s,2)-sin(xi)*'diff(chi,s)+(chi-1)^2 = 
0,'diff(xi,s)+cos(chi)*xi = 0];
ics:['at('diff(chi,s),s = 0)+'at(chi,s = 0) = 0,'at('diff(chi,s),s = 
0)-'at(chi,s = 0) = zeta,'at(xi,s = 0) = 1+zeta];
convert_to_first_order_system(des,[chi,xi],s,ics);