1eval_when([translate,batch,demo,load,loadfile],
2dv(a)::=buildq([a],define_variable(a,'a,any)),
3dvl(a)::=buildq([a],define_variable(a,[],list)),ttyoff:true)$
4
5/*  (c) Copyright 1981 Massachusetts Institute of Technology  */
6
7/* dynamalloc:true$ */
8define_variable(ieqnprint,true,boolean)$
9define_variable(techlist,'[first,all,vlfrnk,transform,collocate,
10		flfrnk2nd,tailor,fredseries,neumann,
11		flfrnk1st,abel,firstkindseries],list)$
12dv(tech)$
13dv(px)$
14dv(eqn)$
15dv(fx)$
16dv(uvar)$
17dv(xvar)$
18dv(kxu)$
19dv(ax)$
20dv(bx)$
21dv(iesoln)$
22dv(napprox)$
23dv(pxlist)$
24dv(eqnlist)$
25dv(inteqn)$
26dv(lowlim)$
27dv(pofx)$
28dv(iclist)$
29dvl(unklist)$
30dv(highlim)$
31dv(%c)$
32dv(eqno)$
33dvl(ueqnlist)$
34dv(x)$
35dv(xhat)$
36dv(opr)$
37dv(guesslist)$
38dv(guess)$
39
40ieqn([arglist]):=
41   if length(arglist)<2 then error("ieqn requires at least two arguments") else
42	ieqn1(arglist[1], arglist[2],
43		if length(arglist)>2 then arglist[3] else [],
44		if length(arglist)>3 then arglist[4] else [],
45		if length(arglist)>4 then arglist[5] else [])$
46
47opr(e):=if mapatom(e) then 'none else inpart(e,0)$
48arg(e):=inpart(e,1)$
49alias(exitfor, return)$
50
51/* identifies functions unknown to macsyma */
52unfun(e):=is(?getchar(e,1)='?\$)$
53
54/* ***** main program *****/
55ieqn1(eqnlist,pxlist,tech,napprox,guesslist):=
56   block([iesoln, xvar, uvar, ax, bx, kxu, fx, eqn, px, guess, system,
57	  inflag, dispflag, singsolve, solveradcan, keepfloat],
58	inflag:true, dispflag:false, system:false,
59	if not listp(eqnlist) then eqnlist:[eqnlist],
60	if not listp(pxlist) then pxlist:[pxlist],
61	if length(eqnlist)#length(pxlist) then
62	  error("number of equations # number of unknowns"),
63	for unk in pxlist do
64	  if mapatom(unk) or length(unk)#1 or not atom(xvar:arg(unk)) or
65	     numberp(xvar) then error(unk," improperly specified unknown"),
66	if tech=[] then myprint("default 3rd arg, technique: ",tech:'first) else
67	if not member(tech,techlist) then error(tech," invalid technique - see techlist"),
68	if napprox=[] then myprint("default 4th arg, number of iterations or coll. parms.: ",napprox:1) else
69	if not integerp(napprox) or napprox<=0 then error("napprox must be a pos. integer"),
70	if guesslist=[] then myprint("default 5th arg, initial guess: ",guesslist:'none),
71	if guesslist='none then guess:'none else
72	  (if not listp(guesslist) then guesslist:[guess:guesslist],
73	   if length(guesslist)#length(pxlist) then error("number of guesses # number of unknowns")),
74	singsolve:solveradcan:keepfloat:true, iesoln:[],
75	eqnlist:maplist(lambda([eqn],num(ratsimp(lhs(eqn)-rhs(eqn)))),eqnlist),
76	if length(pxlist)=1 then (px:first(pxlist), eqn:first(eqnlist))
77	  else system:true,
78	if not system then
79	   (try('vlfrnk), if tech='done then return(reverse(iesoln))),
80	try('transform),
81	if tech='done then return(reverse(iesoln)),
82	for eqnlist in solvesys(eqnlist,pxlist) do
83	   (try('flfrnk2nd), try('tailor),
84	   if not system then try('fredseries), try('neumann)),
85	if not system and firstkind() then
86	   (try('flfrnk1st), try('abel), try('firstkindseries)),
87	try('collocate), return(reverse(iesoln)))$
88
89/* invokes solution method catching and containing errors.  globals: tech */
90try(routine):=if member(tech,['all, 'first, routine]) then
91	block([attempt], attempt:errcatch(catch(apply(routine,[]))),
92	  if attempt=[] then (?princ('? in\ ), ?princ(?stripdollar(routine)),
93		?terpri()) else
94	    if first(attempt) and tech='first then tech:'done)$
95
96
97/* test for 1st kind and extracts parts.
98  globals: eqn, kxu, uvar, xvar, ax, bx */
99firstkind():=block([integral, xpoints, consistent],
100	if not freeof(px,eqn) then return(false),
101	integral:catch(findint(eqn)),
102	if integral=false then return(false),
103	fx:mysolve(eqn,integral), if fx=[] then return(false),
104	fx:first(fx), uvar:inpart(integral,2),
105	kxu:lin(arg(integral), subst(uvar,xvar,px)),
106	if kxu=false then return(false),
107	ax:inpart(integral,3), bx:inpart(integral,4),
108	if kxu[2]#0 then fx:ratsimp(fx-myint(kxu[2],uvar,ax,bx)),
109	xpoints:mysolve(ax-bx,xvar), kxu:kxu[1],
110	if xpoints=[] then return(true) else consistent:false,
111	for xv in xpoints do
112	  if ratsubst(xv,xvar,fx)=0 then exitfor(consistent:true),
113	return(consistent))$
114
115/* returns first integral found in expr */
116findint(expr):=if not mapatom(expr) then
117		(if opr(expr)#'?%integrate then
118			(for i in expr do findint(i), false)
119		 else throw(expr))$
120
121/* globals: ieqnprint, iesoln */
122printsoln(soln,tech,order,kind):=block([dispflag, den],
123	dispflag:ieqnprint, if kind#[] then kind:[kind],
124	if order#[] then kind:cons(order,kind),
125	if not listp(soln) then soln:[soln],
126	soln:maplist(lambda([elmt], block([elem,den],
127	elem:ratsimp(elmt), den:denom(elem),
128	if numberp(den) and den#1 and opr(num(elem))="+" then
129	  multthru(elem) else elem)),soln),
130	if rest(soln)=[] then soln:first(soln),
131	iesoln:cons(?displine(cons(soln,cons(tech,kind))),iesoln),
132	return(true))$
133
134alias(myint, integrate)$
135
136myprint(msg1,msg2):=if ieqnprint then print(msg1,msg2)$
137
138/* if expr = a*var+b returns [a,b] else false */
139lin(expr,var):= block([lc],
140	lc:ratsimp(diff(expr,var)),
141	if freeof(var,lc) then [lc, ratsubst(0,var,expr)])$
142
143/* solves eqn for unk and returns a list of actual solution values */
144mysolve(eqn,unk):=block([result], result:[],
145	eqn:solve(eqn,unk), eqn:maplist('rhs, apply('ev,[eqn])),
146	for ans in eqn do
147	  if freeof(unk, ans) then result:cons(ans, result),
148	return(result))$
149
150/* similar to mysolve but for systems. */
151solvesys(eqns,unks):=
152	(eqns:solve(eqns,unks), eqns:apply('ev,[eqns]), maplist(lambda([el],
153	 if listp(el) then maplist('rhs,el) else [rhs(el)]), eqns))$
154
155/* if expr can be expressed as sum(f[i](x)*g[i](u),i,1,n) then
156  result is [[f1,g1],...[fn,gn]] else [].  globals: xvar  */
157sumfactors(expr,uvar):=block([other, rem],
158	expr:catch(frnk(expr)), if expr=false then return([]),
159	if freeof(uvar,expr) then return([[expr,1]]),
160	if freeof(xvar,expr) then return([[1,expr]]),
161	expr:expand(expr),
162	other:mypartition(expr,uvar), expr:mypartition(expr,xvar),
163	if other[1]<expr[1] then expr:other,
164	rem:if expr[3]=0 then [] else [[expr[3],1]],
165	if expr[4]#0 then rem:cons([1,expr[4]],rem),
166	if expr[2]#0 then for fc in expr[2] do rem:cons(partition(fc,uvar),rem),
167	return(rembox(rem)))$
168
169/* converts e to finiterank form if possible otherwise throws false.
170  globals: xvar, uvar  */
171frnk(e):=
172	if freeof(xvar,e) or freeof(uvar,e) then e else
173	block([op,ar,up], op:opr(e), ar:arg(e),
174	  if member(op,["+","*"]) then return(map('frnk,e)),
175	  if op="^" then (up:inpart(e,2),
176	   if integerp(up) then
177		if up>0 then return(frnk(ar)^up) else throw(false),
178	   up:plussplit(expand(up)),
179	   if up#[] and freeof(uvar,xvar,ar) then
180	       return(box(ar^up[1])*box(ar^up[2])) else throw(false)),
181	  if member(opr(ar),["*","+"]) then (e:partition(ar,xvar),
182	    if not freeof(uvar,e[2]) then throw(false)),
183	  if op='?%log and opr(ar)="*" then log(e[1])+log(e[2]) else
184	  if opr(ar)="+" then
185	    if op='?%sin  then sin(e[1])*cos(e[2])+cos(e[1])*sin(e[2]) else
186	    if op='?%cos  then cos(e[1])*cos(e[2])-sin(e[1])*sin(e[2]) else
187	    if op='?%sinh then sinh(e[1])*cosh(e[2])+cosh(e[1])*sinh(e[2]) else
188	    if op='?%cosh then cosh(e[1])*cosh(e[2])+sinh(e[1])*sinh(e[2]) else
189	    throw(false) else
190	  throw(false))$
191
192/* if e=f(x)+g(u) returns [f(x),g(u)] else [].  globals: xvar, uvar */
193plussplit(e):=
194	if opr(e)="+" then (e:partition(e,uvar),
195	    if freeof(xvar,e[2]) then e else []) else
196	if freeof(uvar,e) then [e,0] else
197	if freeof(xvar,e) then [0,e] else []$
198
199/* solves exs for unks and substitutes in form */
200solveandsubst(exs,unks,form,tech):=block([soln,trial],
201	if (soln:solve(exs,unks))=[] then
202	  (print("for a ",tech," solution substitute in the expression:"),
203	   ldisp(form), apply('print,cons("the values of ",unks)),
204	   print("that make the following expressions simultaneously zero"),
205	   apply('ldisp,exs), return(false)),
206	for sol in apply('ev,[soln]) do (trial:apply('ev,[form,sol]),
207	  apply('printsoln,append([trial,tech],
208	    if tech='collocate then [napprox,testsoln(trial)] else [[],[]]))),
209	return(true))$
210
211/* tests solution for exactness.  globals: eqnlist, pxlist */
212testsoln(resultlist):=block([flag], apply('local,maplist('opr,pxlist)),
213	flag:[], maplist('define,pxlist,resultlist),
214	resultlist:apply('ev,[eqnlist,'integrate,'ratsimp]),
215	for val in resultlist do
216	  if val#0 then exitfor(flag:'approximate),
217	return(flag))$
218
219
220/* fout=h(x,u)+q(x)+r(u).  returns [n,h(x,u),q(x),r(u)] where
221  n is the rank of fout.  globals: xvar, uvar */
222mypartition(fout,var):=block([qx, ru, rem, con],
223	if opr(fout)#"+" or opr(fout:factorout(fout,var))#"+" then
224		return([1,[fout],0,0]),
225	rem:con:qx:ru:0,
226	for fc in fout do
227	  if freeof(uvar,fc) then
228		(if freeof(xvar,fc) then con:con+fc else qx:qx+fc) else
229	  if freeof(xvar,fc) then ru:ru+fc else rem:rem+fc,
230	fout:if rem=0 then 0 else
231	     if opr(rem)="+" then length(rem) else (rem:[rem], 1),
232	if qx#0 then (fout:fout+1, qx:qx+con, con:0),
233	if ru#0 then (fout:fout+1, ru:ru+con),
234	return([fout, rem, qx, ru]))$
235
236/* replaces all integrals and unknown functions in ex with zero  */
237zeroint(ex):=
238	if mapatom(ex) then ex else
239	if opr(ex)='?%integrate or unfun(opr(ex)) then 0 else
240	map('zeroint,ex)$
241
242/* ----------- the following functions apply to 1st or 2nd kind eqns */
243
244/* variable-limit finite-rank 1st or 2nd kind. reduction to ode.
245  globals: eqn, xvar, px */
246vlfrnk():=block([initial, lowlim, unklist, inteqn, pofx, kind,
247			firstkind, difflist, iclist],
248	pofx:px, initial:true, unklist:difflist:[], kind:'incomplete,
249	firstkind:is(freeof(px,eqn)), inteqn:vconvert(eqn),
250	iclist:[xvar=lowlim],
251	for count thru length(unklist) do
252	  ( if firstkind then firstkind:false else
253	      (if initial then getics(), pofx:diff(pofx,xvar)),
254	    difflist:cons(inteqn,difflist), inteqn:diff(inteqn,xvar),
255	    if freeof('?%integrate, inteqn) then exitfor(false)),
256	apply('remove,[opr(px), 'atvalue]),
257	if not freeof('?%integrate, inteqn) then
258	  ( difflist:solve(difflist, unklist),
259	    if difflist=[] then throw(false),
260 inteqn:apply('ev,[inteqn,difflist])),
261	lowlim:derivdegree(inteqn,px,xvar), iclist:reverse(iclist),
262	if lowlim=0 then (iclist:[], lowlim:3,
263	  if (pofx:mysolve(inteqn,px))#[] then (inteqn:first(pofx), kind:[])),
264	if lowlim>2 or (pofx:ode2(inteqn,px,xvar))=false then
265	    return(printsoln(inteqn, 'vlfrnk, iclist, kind)),
266	pofx:ratsimp(pofx),
267	if initial and length(iclist)-1=lowlim and (inteqn:errcatch(
268	   apply(if lowlim=1 then 'ic1 else 'ic2, cons(pofx,iclist))))#[] then
269	     (pofx:first(inteqn), iclist:[]),
270	if opr(pofx)="=" and lhs(pofx)=px then (pofx:rhs(pofx), kind:[]),
271	if iclist=[] and kind#[] then iclist:0,
272	printsoln(pofx, 'vlfrnk, iclist, kind))$
273
274/* obtains initial conditions.
275   globals: inteqn, xvar, lowlim, pofx, iclist, initial */
276getics():=block([val, init],
277	val:at(inteqn, xvar=lowlim), init:at(pofx, xvar=lowlim),
278	init:mysolve(val, init), if init#[] then
279	  ( init:first(init), atvalue(pofx, xvar=lowlim, init),
280	    iclist:cons(pofx=init, iclist)) else initial:false)$
281
282
283/* converts integrands in fun to finiterank form. upper limit must
284be xvar and lower limit must be a constant.
285  globals: initial, lowlim, unklist, xvar */
286vconvert(fun):=
287	if mapatom(fun) then fun else
288	if opr(fun)#'?%integrate then map('vconvert, fun) else
289	if not freeof(xvar,inpart(fun,3)) or
290	   inpart(fun,4)#xvar then throw(false) else
291	block([intgr, newfun, int],
292		if lowlim='lowlim then lowlim:inpart(fun,3) else
293		if lowlim#inpart(fun,3) then initial:false,
294		intgr:sumfactors(arg(fun), inpart(fun,2)),
295		if intgr=[] then throw(false), newfun:0,
296		for term in intgr do
297		  (int:substinpart(term[2], fun, 1),
298		   if not member(int,unklist) then unklist:cons(int, unklist),
299		   newfun:newfun+term[1]*int),
300		return(newfun))$
301
302
303/* laplace transform.   globals: eqnlist, pxlist, xvar  */
304transform():=block([teqnlist, translist, %s, ps, flag],
305	ps:lambda([fun], laplace(fun, xvar, %s)), flag:false,
306	teqnlist:maplist('ps,eqnlist), translist:maplist('ps,pxlist),
307	for soln in solvesys(teqnlist,translist) do
308	 if freeof('?%integrate, '?%laplace, soln) then
309	   (soln:maplist(lambda([fun], ilt(fun,%s,xvar)), soln),
310	   ps:freeof('?%ilt, soln),
311	   printsoln(soln, 'transform, if ps then [] else 0, if ps then [] else 'incomplete),
312	   flag:flag or ps),
313	return(flag))$
314
315/* collocation.  globals: eqnlist, pxlist, napprox, xvar */
316collocate():=block([lowlim, highlim, unklist, %c, elist, form, incr, point,
317			name, listeqns, solnlist],
318	apply('local, cons('%c,maplist('opr,pxlist))), lowlim:'minf,
319highlim:'inf,
320	map('getlimits, eqnlist), if highlim<=lowlim then throw(false),
321	solnlist:listeqns:unklist:[],
322	for unk in pxlist do
323	( name:opr(unk), form:approx(name,napprox,xvar),
324	  apply('define,[unk,form]), solnlist:cons(form,solnlist),
325	  for parm:0 thru napprox-1 do unklist:cons(%c[name,parm],unklist)),
326	elist:apply('ev,[eqnlist,'integrate,'expand]),
327	if not freeof('?%integrate,elist) then throw(false),
328	if freeof(xvar,elist) then listeqns:elist else (point:lowlim,
329	 incr:if napprox>1 then (highlim-lowlim)/(napprox-1) else 0,
330	 for i thru napprox do
331	   (listeqns:append(subst(point,xvar,elist),listeqns),
332	   point:point+incr)),
333	solveandsubst(listeqns,unklist,reverse(solnlist),'collocate))$
334
335/* obtains largest lower limit and least upper limit for collocation points
336	globals: lowlim, highlim */
337getlimits(expr):=
338	if not mapatom(expr) then
339	  if opr(expr)#'?%integrate then
340		for sub in expr do getlimits(sub) else
341	  block([low,high], low:inpart(expr,3), high:inpart(expr,4),
342		if not numberp(low) or not numberp(high) then throw(false),
343		lowlim:max(low,lowlim), highlim:min(high,highlim))$
344
345/* approximation function for unknown function fun(var)  */
346approx(fun,nparms,var):=sum(%c[fun,i-1]*var^(i-1),i,1,nparms)$
347
348/* ---- the following functions apply only to 2nd kind eqns */
349
350/*  fixed-limit, finite-rank, 2nd kind  globals:  eqnlist, pxlist  */
351flfrnk2nd():=block([frlist, unklist, ueqnlist, eqno, %c],
352	apply('local, cons('%c,maplist('opr, pxlist))),
353	unklist:ueqnlist:[], eqno:0,
354	frlist:maplist('fconvert, eqnlist),
355	maplist('define, pxlist, frlist),
356	ueqnlist:apply('ev,[ueqnlist,'integrate]),
357	solveandsubst(ueqnlist, unklist, frlist, 'flfrnk2nd))$
358
359/*  returns result of replacing all occurences of
360	'integrate(f(x,u),u,a,b) in fun with
361	sum(q[j](x)*%c[j],j,1,n) where %c[j] is 'integrate(r[j](u),u,a,b)
362	after expressing integrands in finite-rank form
363	(here f(x,u) becomes sum(q[j](x)*r[j](u),j,1,n).
364	globals: eqno (current number of subscripted unknowns %c),
365		 ueqnlist (list of equations relating %c's)
366		 unklist (list of all %c unknowns to be solved for)
367		 xvar (indepedent variable) */
368fconvert(fun):=
369	if mapatom(fun) then fun else
370	if opr(fun)#'?%integrate then map('fconvert, fun) else
371	if not freeof(xvar,inpart(fun,3)) or not freeof(xvar, inpart(fun,4))
372	 then throw(false) else
373	block([newfun, intgrnd],
374	  intgrnd:sumfactors(arg(fun), inpart(fun,2)),
375	  if intgrnd=[] then throw(false),  newfun:0,
376	  for term in intgrnd do
377		(eqno:eqno+1, newfun:newfun+%c[eqno]*term[1],
378		 ueqnlist:cons(%c[eqno]-substinpart(term[2],fun,1), ueqnlist),
379		 unklist:cons(%c[eqno], unklist)),
380	  return(newfun))$
381
382/*  taylor series. globals: eqnlist, pxlist, xvar, napprox  */
383tailor():=block([xhat, ufun, eqtn, tfun, neqns, vlist, order, value, fact],
384	apply('local,append([ufun,eqtn,tfun],maplist('opr,pxlist))),
385	map('getxhat, eqnlist), neqns:0, vlist:pxlist,
386	for expr in eqnlist do
387	( neqns:neqns+1, eqtn[neqns]:expr, ufun[neqns]:first(vlist),
388	  vlist:rest(vlist), tfun[neqns]:value:subst(xhat,xvar,expr),
389	  atvalue(ufun[neqns], xvar=xhat, value)),
390	fact:1, order:napprox,
391	for i thru napprox do
392	( fact:fact*i,
393	  for j thru neqns do
394	  ( value:errcatch(diff(eqtn[j],x)),
395	    if value=[] then exitfor(order:i-1),
396	    eqtn[j]:first(value), value:errcatch(ratsimp(at(eqtn[j],xvar=xhat))),
397	    if value=[] then exitfor(order:i-1),
398	    ufun[j]:diff(ufun[j],x), value:first(value),
399	    atvalue(ufun[j], xvar=xhat, value),
400	    tfun[j]:tfun[j]+value*(x-xhat)^i/fact),
401	  if order#napprox then exitfor(false)),
402	for i:neqns step -1 thru 1 do vlist:cons(tfun[i],vlist),
403	printsoln(vlist, 'tailor, order, testsoln(vlist)))$
404
405/* obtains expansion point for taylor series by solving b(x)=a(x).
406  globals: xvar, xhat */
407getxhat(expr):=
408	if mapatom(expr) then false else
409	if opr(expr)#'?%integrate then for sub in expr do getxhat(sub) else
410	if xhat='xhat then
411	  (xhat:mysolve(inpart(expr,3)-inpart(expr,4),xvar),
412	   if xhat=[] then throw(false) else xhat:first(xhat),
413	   if ratsubst(xhat, xvar, inpart(expr,3))#xhat then throw(false)) else
414	if subst(xhat,xvar,expr)#0 then throw(false)$
415
416/* fredholm-carleman series.   globals: napprox, kxu, ax, bx, fx, xvar, uvar */
417fredseries():=block([order, alpha, bta, top, bot, kind, tvar, eqtn, kxt],
418	eqtn:first(eqnlist),
419	alpha:catch(findint(eqtn)), if alpha=false then return(false),
420	uvar:inpart(alpha,2), bta:lin(eqtn,alpha),
421	if bta=false then return(false),
422	kxu:lin(arg(alpha),subst(uvar,xvar,px)),
423	if kxu=false then return(false),
424	ax:inpart(alpha,3), bx:inpart(alpha,4), fx:bta[2],
425	if kxu[2]#0 then fx:fx+myint(bta[1]*kxu[2],uvar,ax,bx),
426	kxu:kxu[1]*bta[1], order:napprox, kind:'approximate,
427	kxt:subst(tvar,uvar,kxu), bot:bta:1,
428	top:alpha:rat(kxu,uvar,xvar),
429	for i thru napprox do
430	( bta:-myint(ratsubst(uvar,xvar,alpha),uvar,ax,bx)/i,
431	  alpha:bta*kxu+myint(kxt*ratsubst(tvar,xvar,alpha),tvar,ax,bx),
432	  bot:bot+bta,
433	  if alpha=0 then exitfor(kind:[]),
434	  top:top+alpha),
435	kxt:fx+myint(subst(uvar,xvar,fx)*top,uvar,ax,bx)/ratdisrep(bot),
436	printsoln(kxt, 'fredseries, order, kind))$
437
438
439/* neumann series.   globals: eqnlist, pxlist, guesslist, napprox  */
440neumann():=block([order, newguess],
441	apply('local, maplist('opr,pxlist)), order:napprox,
442	if guesslist='none then guesslist:maplist('zeroint, eqnlist),
443	for count thru napprox do
444	  (maplist('define, pxlist, guesslist),
445	   newguess:apply('ev,[eqnlist,'integrate]),
446	   if not freeof('?%integrate, newguess) then exitfor(order:count-1),
447	   guesslist:maplist('ratsimp,newguess)),
448	printsoln(guesslist, 'neumann, order, testsoln(guesslist)))$
449
450
451
452/* ---- the following functions apply only to 1st kind eqns */
453
454/* fixed-limit finite-rank 1st kind.  globals: kxu, fx, uvar, xvar, ax, bx */
455flfrnk1st():=block([sf, cnt, intg, form, %c, unklist, eqlist,res],
456
457	local(%c), cnt:form:0, unklist:[],
458	if not freeof(xvar,[ax,bx]) then return(false),
459	intg:sumfactors(kxu,uvar),
460	for term in intg do
461	  (cnt:cnt+1, form:form+term[2]*%c[cnt],
462	   unklist:cons(%c[cnt],unklist)),
463	eqlist:[],
464	sf:num(ratsimp(myint(kxu*form,uvar,ax,bx)-fx,xvar)),
465	res:0, if opr(sf)#"+" then sf:[sf],
466	for term in sf do
467	  if opr(term)="*" then
468	    (intg:partition(term,xvar),
469	     if intg[2]=1 then res:res+term else eqlist:cons(intg[1],eqlist))
470	  else if freeof(xvar,term) then res:res+term else error("error 1"),
471	if res#0 then eqlist:cons(res,eqlist),
472	if length(eqlist)#cnt then error("error 2"),
473	solveandsubst(eqlist, unklist, subst(xvar,uvar,form), 'flfrnk1st))$
474
475/* some=one*peace+two, where peace could be a product. */
476mydivide(some,peace,var):=block([res, one, two],
477	one:two:0, if opr(some)#"+" then some:[some],
478	for prt in some do
479	  (res:quotient(prt,peace),
480	   if res#0 and freeof(var,res) then one:one+res else two:two+prt),
481	return([one,two]))$
482
483
484/* generalized abel method.  globals: kxu, ax, bx, fx, xvar, uvar */
485abel():=block([power,den,fun],
486	if not freeof(xvar,ax) or bx#xvar or opr(kxu)#"^" then throw(false),
487	power:-inpart(kxu,2), den:inpart(kxu,1),
488	if not freeof(uvar, power) or opr(den)#"+" then throw(false),
489	fun:partition(den,uvar),
490	if not freeof(xvar, fun[2]) or
491	   ratsimp(subst(uvar,xvar,fun[1])+fun[2])#0 then throw(false),
492	if sign(power)#'pos or sign(power-1)#'neg then throw(false),
493	fun:myint(diff(-fun[2],uvar)*subst(uvar,xvar,fx)*den^(power-1),uvar,ax,bx),
494	fun:sin(%pi*power)/%pi*diff(fun,xvar),
495	printsoln(fun, 'abel, [], []))$
496
497
498/* firstkindseries.  globals: px, napprox, guess, eqn, px  */
499firstkindseries():=block([kind, order, correction, trial],
500	apply('local,[opr(px)]), kind:'approximate, order:napprox,
501	trial:if guess='none then zeroint(eqn) else guess,
502	for i thru napprox do
503	  (apply('define,[px, trial]), correction:apply('ev,[eqn,'integrate]),
504	   if not freeof('?%integrate, correction) then exitfor(order:i-1),
505	   if (correction:ratsimp(correction))=0 then exitfor(kind:[]),
506	   trial:ratsimp(trial-correction)),
507	printsoln(trial, 'firstkindseries, order, kind))$
508
509/* kill(labels)$ */
510ttyoff:false$
511