#==============================================================================
# Utility procedures for working with expressions that contain derivatives.
#==============================================================================
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washington, DC 20392
# murison@aa.usno.navy.mil
# http://aa.usno.navy.mil/murison/
#==============================================================================

with(utils):

ode[init] := proc()
	global `type/diff`;
	`type/diff` := proc(x)
		isdiff(x);
	end proc;
end proc;


#==============================================================================
# Determine if an expression is a derivative.
#==============================================================================
ode[isdiff] := proc(x)
	if type(x,function) and op(0,x)='diff' then
		true;
	else
		false;
	end if;
end proc;



#==============================================================================
# Determine the order of a derivative.
#==============================================================================
ode[difforder] := proc( diffexpr::diff )

	local n, d;

    n := 0;
    d := diffexpr;
    while type(d,'diff') do
        n := n+1;
        d := op(1,d);
    od;

    return n;
end proc;



#==============================================================================
# Return the order of the highest-order
# derivative in an expression or equation.
#==============================================================================
ode[maxPDEorder] := proc( expr::{algebraic,`=`,indexed,list,array} )

    local p, n;

    if type(expr,`=`) then
        return  procname( lhs(expr)-rhs(expr) );
    end if;

	if isdiff(expr) then
		return  difforder(expr);
	end if;

	if nops(expr)=1 then
		return  0;
	end if;

    n := 0;
    for p in expr do
        if isdiff(p) then
            n := max( n, difforder(p) );
		else
			n := max( n, procname(p) );
        end if;
    od;

    return n;
end proc;



#==============================================================================
# Return the variables of a derivative.  The return
# quantity is a list of the form [var0,var1,...,varn]
# where var0 is the dependent variable and var1..varn are
# the n independent variables, where n = order of diffexpr.
# Hence, diff( y(x,t), x$2, t ) would produce z := [y(x,t),x,x,t].
# Then diff( op(z) ) would recover the original derivative.
#==============================================================================
ode[diffvars] := proc( diffexpr::diff )

	local i, dvar, ivars;

	ivars := [];
	dvar  := diffexpr;
	for i from 1 to difforder(dvar) do
		ivars := [op(ivars),op(2,dvar)];
		dvar  := op(1,dvar);
	od;

	return  [ dvar, op(ivars) ];
end proc;



#==============================================================================
# Return the dependent variable of a derivative.
#==============================================================================
ode[depvar] := proc( diffexpr::diff )
	local vars;
	vars := diffvars( diffexpr );
	return  vars[1];
end proc;



#==============================================================================
# Return the coefficients of the derivatives in the expression expr,
# which can be an algebraic expression or an equation.  The returned
# structure is of the form L = [ [c1,c2,...,cn], [d1,d2,...,dn] ], where
# L[1] is the list of coefficients and L[2] is the list of corresponding
# derivatives.  Hence, the original ODE can be recovered with
# sum( L[1][k]*L[2][k], k=1..nops(L[1]) ) = 0;
#==============================================================================
ode[diffcoeffs] := proc( expr::{algebraic,algebraic=algebraic} )

	local ODE, p, dlist, clist, parts, i, L, L0;

	if type(expr,algebraic=algebraic) then
		ODE := lhs(expr) - rhs(expr);
	elif type(expr,algebraic) then
		ODE := expr;
	end if;

	if not has(ODE,diff) then
		ERROR(`: input contains no derivatives`);
	end if;

	#
	# Collect on diffs first.
	#
	ODE := collect(ODE,diff,distributed);

	#
	# Special case: form coeff*diff
	#
	if type(ODE,`*`) then
		parts := {op(ODE)};
		dlist := select( has, parts, diff );
		clist := remove( has, parts, diff );
		return  [ [convert(clist,`*`)], [convert(dlist,`*`)] ];
	end if;

	#
	# Special case: form (diff)^n, n = [1...]
	#
	if type(ODE,`^`) or isdiff(ODE) then
		return  [[1],[ODE]];
	end if;

	#
	# General case.  Separate coefficients and corresponding derivatives.
	#
	dlist := [];	#list of derivatives
	clist := [];	#corresponding list of coefficients
	parts := [op(ODE)];
	L0    := [];
	for p in parts do
		if has(p,'diff') then
			L     := procname(p);
			clist := [op(clist),op(L[1])];
			dlist := [op(dlist),op(L[2])];
		else
			L0 := [ op(L0), p ];
		end if;
	od;
	if nops(L0) > 0 then
		clist := [op(clist),convert(L0,`+`)];
		dlist := [op(dlist),1];
	end if;

	return  [ clist, dlist ];
end proc;



#==============================================================================
# Replace the coefficients of derivatives in an equation
# with the dummy coefficients A[k], where k={1..N} and N is
# the total number of coefficients (i.e., the number of derivative terms).
# The output is a table T of the form
#	T['eqn']		The equation(s) with coefficient substitutions made
#	T['coeffs']		A list of the actual coefficients, ordered corresponding
#					to the A[k]
#	T['derivs']		A list of the derivatives corresponding to the A[k]
#   T['vars']		A list of derivative variables [see diffvars()]
# The original equation may be reconstructed by feeding the output of
# MakeCoeffs() into the routine Reconstruct().
#
# Side effects: The dummy coefficient table A[k] is a global variable.
# Anything you might have assigned to the global variable 'A' will be lost.
#==============================================================================
ode[MakeCoeffs] := proc( expr::{algebraic,algebraic=algebraic} )

	local	m, L, c, d, ans, v, p, tmp, q, M;
	global	A;

	L := diffcoeffs(expr);
	c := L[1];				#coefficients
	d := L[2];				#derivatives
	v := [];				#variables
	for p in d do
		if type(p,`^`) then
			M := procname(op(1,p));
			v := [op(v),op(M['vars'])];
		elif type(p,`*`) then
			tmp := [];
			for q in p do
				M   := procname(q);
				tmp := [op(tmp),op(1,M['vars'])];
			od;
			v := [op(v),tmp];
		elif p=1 then
			v := [op(v),[]];
		elif isdiff(p) then
			v := [op(v),diffvars(p)];
		else
			ERROR(`: can't figure out what p is!`);
		end if;
	od;

	A := 'A';
	ans['eqn']    := sum( A[m]*d[m], m=1..nops(c) );
	ans['coeffs'] := c;
	ans['derivs'] := d;
	ans['vars']   := v;

	return  eval(ans);
end proc;



#==============================================================================
# Get the location k of coefficient A[k] in expr, where
# expr is given by MakeCoeffs(), for the derivative diff( var, indvar$ord ).
#==============================================================================
ode[GetCoeffLoc] := proc( expr, var, ord, indvar )

	local p, q, loc;

	# find the term(s) with the desired derivative
	select( has, expr[eqn], diff(var,indvar$ord) );

	if type(%,`+`) then	#more than one term found
		convert(%,set);
	end if;

	select( has, %, A );

	if nops(%)>1 then	#selection is a set of more than one
		q := %;
		# return a set of coefficients
		for p in q do
			loc := location( q, p );
			p   := select(has,p,A);
			q   := subsop( loc=op(1,p), q );
		od;
		q;
	else
		op(1,%);
	end if;
end proc;



#==============================================================================
# Get the coefficient A[k] in expr, where
# expr is given by MakeCoeffs(), for the derivative diff( var, indvar$ord ).
#==============================================================================
ode[GetCoeff] := proc( expr, var, ord, indvar )
  expr[coeffs][GetCoeffLoc(expr,var,ord,indvar)];
end proc;



#==============================================================================
# Substitute the coefficients A[k] as given by MakeCoeffs() for the actual
# coefficients in an equation or equations containing derivatives.  The
# input expr must be a table as output by MakeCoeffs().  See MakeCoeffs()
# for further information.
#==============================================================================
ode[Reconstruct] := proc( expr )

	local k, n, f;

	f := copy(expr['eqn']);
	n := 1;
	for k from 1 to nops(expr['coeffs']) do
		f := subs( A[n]=expr['coeffs'][n], f );
		n := n + 1;
	od;

	return  eval(f);

end proc;



#==============================================================================
# Collect on derivatives.  Optional second argument is a function to 
# apply to the collected terms.  Default operation on collected terms 
# is 'factor'.  To specify no function, use the argument 'nothing'.
#==============================================================================
ode[CollectDiffs] := proc( ODE::{`=`,algebraic,set,list,array} )

	local M, ord, pwr, tlist, slist, k, p, q,
		  expr, n, S, fargs, _H_, _A_;

	if type( ODE, {`=`,set,list,array} ) then
		return  map( procname, args );
	end if;

	if not hastype( ODE, diff ) then
		return ODE;
	end if;

	M   := MakeCoeffs(ODE);

	#
	# sort the equation by order and power
	#
	if not type( M['eqn'], `+` ) then
		slist := [M['eqn']];
	else
		slist := sort( convert(M['eqn'],list), `CollectDiffs/sort` );
	end if;
	ord   := -1;
	pwr   := -1;
	tlist := [];
	k     := 1;
	for p in slist do
		q := select(has,p*_A_,diff);
		if q=NULL then
			q:=1;
		end if;
		# order/power change -- time for a new expression group
		if ord<>maxPDEorder(q) or pwr<>degree(q) then
			ord     := maxPDEorder(q);
			pwr     := degree(q);
			expr[k] := convert(tlist,`+`);
			tlist   := [];
			k       := k + 1;
		end if;
		tlist := [op(tlist),p];
	od;
	expr[k] := convert(tlist,`+`);

	n := k;
	S['coeffs'] := M['coeffs'];
	for k from 1 to n do
		S['eqn'] := copy(expr[k]);
		expr[k]  := Reconstruct(S)*_H_^(n+1-k);
	od;

	if nargs > 1 then
		fargs := args[2..nargs];
	else
		fargs := `factor`;
	end if;

	sum(expr['k'],'k'=1..n);
	collect( %, _H_, fargs );
	`CollectDiffs/collect`(%,fargs);
	expr := subs( _H_=1, % );

	if factor(simplify(expr-ODE)) <> 0 then
		debug_print(procname,`FAILED ERROR CHECK!`,0);
		debug_print(procname,`collected expr not equal to input ODE!`,
					0,`expr - ODE`=%);
		ERROR(`: Quitting.`);
	end if;

	expr;

end proc;


ode[`CollectDiffs/collect`] := proc(expr)
	local p, tmp, loc, newexpr, fargs, q;
	fargs   := args[2..nargs];
	newexpr := expr;
	for p in newexpr do
		if has(p,diff) then
			tmp := select(has,p,diff);
			if has(tmp,diff) then
				loc := location(newexpr,tmp);
				if nops(loc)>0 then
					tmp := collect(tmp,diff,fargs);
					newexpr := subsop( loc=tmp, newexpr );
				elif type(tmp,`*`) then
					for q in tmp do
						loc := location(newexpr,q);
						if nops(loc)>0 then
							q := collect(q,diff,fargs);
#							q := procname(q,fargs);
							newexpr := subsop( loc=q, newexpr );
						end if;
					od;
				end if;
			end if;
		end if;
	od;
	newexpr;
end proc;


#ode[`CollectDiffs/collect`] := proc(expr)
#	local p, tmp, loc, newexpr, fargs;
#	fargs   := args[2..nargs];
#	newexpr := expr;
#	for p in newexpr do
#		tmp := select(has,p,diff);
#		if has(tmp,diff) then
#			loc := location(newexpr,tmp);
#			if nops(loc)>0 then
#				tmp := collect(tmp,diff,fargs);
#				newexpr := subsop( loc=tmp, newexpr );
#			end if;
#		end if;
#	od;
#	newexpr;
#end proc;


ode[`CollectDiffs/sort`] := proc(p,q)
	local a, b, A;
	a := select(has,p*A,diff);
	b := select(has,q*A,diff);
	if a=NULL then a:=1 end if;
	if b=NULL then b:=1 end if;
	if maxPDEorder(a) = maxPDEorder(b) then
		if degree(a) > degree(b) then
			true;
		else
			false;
		end if;
	elif  maxPDEorder(a) > maxPDEorder(b) then
		true;
	else
		false;
	end if;
end proc;



#==============================================================================
# Transform coordinates in a set of ODEs according to the
# transform set xform, then manipulate the resulting transformed
# differential equations until the second-order derivatives appear
# as isolated terms.  Returned the transformed, isolated, differential
# equations as a list.  Optional third argument is the order of
# the derivatives to be isolated (default=2).
#==============================================================================
ode[IsolateDerivs] := proc( ODEs::{set,list}, xform::{set,list} )

	local  T, c, d, globeq, k, M, tmp, depvars, i, j, dset,
				 ODEset, p, q, indvar, n, subslist, eq, ord;
	global time0;

	time0 := time();

	if nargs=2 then
		ord := 2;
	elif nargs=3 then
		ord := args[3];
	else
		ERROR(`: wrong number of arguments`);
	end if;

	#
	# transform the ODEs
	#
	debug_print(procname,`transforming the equations...`,0);
	ODEset := varchange( ODEs, xform );
	debug_print(procname,`transformed equations:`,4,ODEset);
	n := nops(ODEset);

	#
	# mult each eqn by a parameter T[k] and add together
	#
	debug_print(procname,`forming global equation...`,0);
	for k from 1 to n do
		debug_print(procname,`simplifying equation `.k.`...`,1);
		tmp[k] := collect(ODEset[k],diff,factor)*T[k];
	od;
	globeq := sum(tmp['k'],'k'=1..n);
	debug_print(procname,`global equation:`,4,globeq);

	#
	# decompose global eqn
	#
	debug_print(procname,`decomposing global equation...`,0);
	M := MakeCoeffs(globeq);
	debug_print(procname,`checking for decomposition errors...`,0);
	Reconstruct(%);
	tmp := factor(expand(%-(lhs(globeq)-rhs(globeq))));
	debug_print(procname,`error check:`,0,tmp);

	#
	# extract dependent and independent vars
	#
	debug_print(procname,`global differential equation form:`,1,M[eqn]);
	depvars := {};
	indvar  := {};
	for p in M[vars] do
		if nops(p)>0 then
			if type(p[nops(p)],list) then
				for q in p do
					depvars := depvars union convert(q,set);
					indvar  := indvar union {q[nops(q)]};
				od;
			else
				depvars := depvars union convert(p,set);
				indvar  := indvar union {p[nops(p)]};
			end if;
		end if;
	od;
	if nops(indvar)<>1 then
		ERROR(`: more than one independent variable!`);
	else
		tmp    := indvar[1];
		indvar := 'indvar';
		indvar := tmp;
	end if;
	depvars := depvars minus {indvar};
	if nops(depvars)<>n then
		ERROR(`: number of dependent vars does not equal number of ODEs!`);
	end if;
	debug_print(procname,`independent variable:`,0,indvar);
	debug_print(procname,`dependent variables:`,0,depvars);

	# get the coefficients of the indicated derivs
	debug_print(procname,`getting order `.ord.` coefficients...`,0);
	for k from 1 to n do
		c[k] := GetCoeff(M,depvars[k],ord,indvar);
	od;
	debug_print( procname,
				`coefficients of order `.ord.` derivatives:`, 3,
				 convert(c,set) );

	#
	# for each depvar, solve for coefficients that will isolate
	# the derivative of order ord in the global differential equation
	#
	debug_print(procname,`solving coefficients to isolate derivatives...`,0);
	for k from 1 to n do

		debug_print( procname,
					`isolating derivative `.k.`...`, 1,
					 diff(depvars[k],indvar$ord) );

		# manipulate copies of the c[]
		for j from 1 to n do
			d[j] := c[j];
		od;
		d[k] := c[k]=1;  #this is the one we'll isolate
		debug_print(procname,`isolating coefficient c[`.k.`]:`,3,c[k]);

		# solve for the T[k]
		debug_print(procname,`solving for the T[k]...`,2);
		dset     := convert(d,set);
		subslist := solve( dset, {seq(T[j],j=1..n)} );
		if subslist={} then
			debug_print( procname,
						`no solutions for the T[k] in`, 0, dset );
			ERROR(`: Quitting from solve for T[k] section.`);
		end if;
		# solve() converts hyperbolic trig to exp, so we have to undo that
		if has(subslist,exp) then
			subslist := map( convert, subslist, trig );
		end if;
		subslist := map( simplify, subslist );
		debug_print(procname,`coefficient substitutions list:`,2,subslist);

		# subs the T[k] solutions into the global equation to
		# isolate the desired derivative of order ord
		debug_print(procname,`simplifying isolated equation `.k.`...`,1);
		collect( subs(subslist,globeq), diff, simplify );
		eq[k] := CollectDiffs(%);
		debug_print(procname,`isolated equation `.k.`:`,4,eq[k]);

	od;

	debug_print(procname,`done!`,0);
	[seq(eq[i],i=1..n)];

end proc;



#==============================================================================
# Transform the dependent variables of a univariate differential equation.
#==============================================================================
ode[varchange] := proc( deqns::{algebraic,`=`,list,set}, trans::{`=`,set,list} )

	local p, eqs, neweqs;

	#
	# error checks
	#
	if type(deqns,{list,set}) then
		for p in deqns do
			if not has(p,diff) then
				ERROR(`: one or more of deqns is not a differential equation or expression`);
			end if;
		od;
	elif not has(deqns,diff) then
		ERROR(`: deqns is not a differential equation or expression`);
	end if;

	#
	# loop over each differential equation
	#
	eqs := deqns;
	if not type(deqns,{list,set}) then
		eqs := [eqs];
	end if;
	neweqs := [];
	for p in eqs do
		p := subs( trans, p );
		neweqs := [ op(neweqs), p ];
	od;

	if nops(neweqs)=1 then
		op(neweqs);
	else
		neweqs;
	end if;
end proc;



#==============================================================================
# Transform independent variable t to new independent variable T in
# the differential equation deqn.  Teqn must be an equation of the form
#   g(T) = f(t)
# where T is the new independent variable, f is some function of t,
# and g is some function of T.
#==============================================================================
ode[indvarchange] := proc( deqn, Teqn::equation, t::name, T::name )

	local	s, k, tsol, eqn, unknowns, p;
	global	_sol_, _Teqn_, _t_;

	#
	# Annoying error checks.
	#
	if not has( deqn, t ) then
		ERROR(`: deqn does not contain '`.t.`' as independent variable`);
	end if;

	if not has( Teqn, t ) then
		ERROR(`: Teqn does not contain '`.t.`'`);
	elif not has( Teqn, T ) then
		ERROR(`: Teqn does not contain '`.T.`'`);
	end if;

	if nargs=5 and not type(args[5],{set,list}) then
		ERROR(`: unknown functions must be a Maple list or set`);
	end if;

	#
	# Solve Teqn for t.
	#
	tsol := `indvarchange/solve`( Teqn, t );

	#
	# Now make the independent variable change.
	#
	unknowns := {NULL};
	if nargs=5 then
		unknowns := args[5];
	end if;
	# need negative copies to catch unary minus terms in deqn
	for p in unknowns do
		unknowns := unknowns union {-p};
	od;
	eqn := `indvarchange/varchange`( deqn, Teqn, t, T, tsol, unknowns );

	# cleanup of indvarchange/solve globals
	_sol_  := '_sol_';
	_Teqn_ := '_Teqn_';
	_t_    := '_t_';

	eqn;

end proc;


ode[`indvarchange/solve`] := proc( Teqn, t )

	local	k, tt, s;
	global	_sol_, _Teqn_, _t_;

	s := solve( Teqn, t );
	if nops([s])=0 then
		ERROR(`: no solution found for `.t.``);
	elif nops([s]) > 1 then
		if nops(_Teqn_) > 1 and lhs(_Teqn_)=lhs(Teqn)
		   and rhs(_Teqn_)=rhs(Teqn) and _t_=t then
			_sol_;
		else
			debug_print(procname,`nonunique solution for `.t.`:`,0,s);
			debug_print( procname,
						`choose a branch, followed by a `.
						`semicolon (e.g., "2;"): `, 0 );
			k := readstat(`==> `);
			if k < 1 or k > nops([s]) then
				ERROR(`: Quitting.`);
			end if;
			print(s[k]);
			_Teqn_ := Teqn;			# save for future identical calls
			_t_    := t;
			_sol_  := s[k];
		end if;
#	elif type(s,RootOf) then
#		debug_print( procname, `Solution for `.t.` yields a RootOf`, 0, s );
#		tt := simplify( s, assume=real );
#		debug_print( procname, `simplify/assume=real yields`, 0, tt );
#		if type(tt,RootOf) then
#			ERROR(`: Quitting.`);
#		end if;
#		debug_print( procname,
#					`choose a solution, followed by a `.
#					`semicolon (e.g., "1;"): `, 0 );
#		k := readstat(`==> `);
#		if k < 1 or k > nops(tt) then
#			ERROR(`: Quitting.`);
#		end if;
#		print(tt[k]);
#		if type(tt[k],RootOf) then
#			ERROR(`: Quitting.`);
#		end if;
#		tt[k];
	else
		s;
	end if;
end proc;


ode[`indvarchange/varchange`] := proc( deqn, Teqn, t, T, tsol, unknowns )

	local p, Z, g, flist, k, parts, diffs, nodiffs, loc;

	Z := deqn;

	#
	# special case: deqn is a diff
	#
	if isdiff(Z) then
		g := depvar(Z);
		Z := `indvarchange/diff`( subs(t=T,g), Teqn, t, T, difforder(Z) );
		Z := `indvarchange/subs/nondiffs`( t=tsol, Z );
		return Z;
	end if;

	#
	# work on the parts of Z, one at a time
	#
	parts := {op(Z)};

	#
	# split the parts up into terms with and without differentials
	#
	diffs   := select( has, parts, diff );
	nodiffs := parts minus diffs;

	# need negative copies of nodiffs to catch unary minus terms in Z
	for p in nodiffs do
		nodiffs := nodiffs union {-p};
	od;

	#
	# make a list, flist, of the "unknown" functions that occur
	# in the non-differential parts
	#
	flist := select( type, nodiffs, function ) intersect unknowns;

	#
	# recursively do the parts
	#
	for k from 1 to nops(parts) do
		p   := parts[k];
		loc := location( Z, p );
		#skip if this term canceled out and no longer is in Z
		if nops(loc) > 0 then
			if isdiff(p) then
				g := depvar(p);
				p := `indvarchange/diff`( subs(t=T,g), Teqn, t, T, difforder(p) );
				p := `indvarchange/subs/nondiffs`( t=tsol, p );
				Z := subsop( loc=p, Z );
			elif p=t then
				Z := subsop( loc=tsol, Z );
			elif has(p,t) then
				if nops(flist) > 0 and (member(p,flist) or member(-p,flist)) then
					Z := subsop( loc=subs(t=T,p), Z );
				else
					p := procname( p, Teqn, t, T, tsol, unknowns );
					Z := subsop( loc=p, Z );
				end if;
			end if;
		end if;
	od;

	return Z;
end proc;



#==============================================================================
# Perform a change of independent variable differentiation.  Given
# some function g(T), we wish to transform the derivative of g wrt
# the "old" independent variable, t, to the derivative of g wrt a
# "new" independent variable, T, where T = f(t).  Hence,
#
#   dg/dt = (df/dt)*dg/dT
#
# Teqn is the transformation of independent variables from T to t and
# is of the form T = f(t), where f is some function of the "old"
# independent variable, t.  The result will in general be an
# expression containing both t and T.  An optional fourth argument is
# the number of times to perform the differentiation.
#==============================================================================
ode[`indvarchange/diff`] := proc( expr, Teqn::equation, t::name, T::name )

	local f, newexpr;

	f := `indvarchange/solve`( Teqn, T );

	if not has(expr,T) then
		ERROR(`: expr does not contain the new independent variable `.T);
	end if;
	if not has(f,t) then
		ERROR(`: rhs of Teqn does not contain old independent variable `.t);
	end if;

	newexpr := expr;
	if nargs=5 and args[5] > 1 then
		newexpr := procname( newexpr, Teqn, t, T, args[5]-1 );
	end if;

	if f = t then	#special case: Teqn := T=t -- treat T as a function of t
		diff(newexpr,t) + diff(T(t),t)*diff(newexpr,T);
	else			#Teqn := T=f(t)
		diff(newexpr,t) + diff(f,t)*diff(newexpr,T);
	end if;
end proc;



#==============================================================================
# Perform subs() on the parts of expr that are not derivatives
#==============================================================================
ode[`indvarchange/subs/nondiffs`] := proc( s, expr )

	local parts, p, loc, newexpr;

	newexpr := expr;

	if not isdiff(expr) then
		parts   := [op(expr)];
		newexpr := expr;
		for p in parts do
			#skip if p has nothing useful to work with
			if nops( indets(p) ) > 0 then
				#skip if p is a derivative
				if not isdiff(p) then
					loc := location( newexpr, p );
					#skip if this term canceled out and no longer is in newexpr
					if nops(loc) > 0 then
						if nops(p)=1 then
							p := subs( s, p );
						else
							p := procname( s, p );
						end if;
						newexpr := subsop( loc=p, newexpr );
					end if;
				end if;
			end if;
		od;
	end if;

	newexpr;
end proc;



save( ode, cat(ODE_PATH,"ode.m") );