/* ************** beginning of main ******************************** */

lode2(oeq):=block([lhs,ay,ypx],tobegin(dum),
   remvalue(y1,y2), swro:0, ys:f, nhstddeq(oeq), if rr#0 then 
    (print("first",str1,eq)), kt:0, ys:hlode2(eq), 
   if rr=0 then return(ys), print("the solution of the homog. eq. is",ys),
   ldf0(8), ypx:nonhom(rr,y1,y2), return(ys+ypx))$

tobegin(dum):=block([],dyx2:'diff(y,x,2),dyx1:'diff(y,x),
   dyt2:'diff(y,t,2),dyt1:'diff(y,t),dux2:'diff(u,x,2),dux1:'diff(u,x),
   str1:"we solve",str2:"we use",str3:"the result is",str4:"type y or n",
   str5:"positive,zero,or negative? type p,z,or n")$

nhstddeq(oeq):=block([], weq:expand(oeq), lhs:part(weq,1), 
   rhs:part(weq,2), lhs:lhs-rhs, deqcoeff(lhs), wlhs:lhs, 
   lhs:za*dyx2+zb*dyx1+zc*y, rhs:-ratsimp(wlhs-lhs), 
   weq:lhs=rhs, steq:stddeq(weq), lhs:part(steq,1), eq:lhs=0,
   if lhs#wlhs then print(str1,steq))$

hlode2(eq):=block([exp,ro], aeq[kt]:eq, swsv:f, 
 ll1, leq:part(eq,1), tp:brchde1(leq),
   if zc=0   then (ldf0( 9),svode1(rb),return(ys)),
   if tp= 50 then (ldf0(10),skelp(rb,rc), if ys#f then return(ys)),
   if swpt=y then (ldf0(11),gptm(rb,rc),  if ys#f then return(ys) 
      else if rptm=t then (if freeof(y,vtr) then idvprep(1) 
                           else idvprep(2),go(ll1))),
   if tp>60  then (lv:listofvars(za), if length(lv)>1 then 
      (ldf0(16), if stddsp(za)=t then go(ll2)), tblspt(za,rb,rc,x,n), 
      if nspt>3  then (ldf0(12),if grmaps(dum)=y then go(ll1))),
   ha:hipow(za*q1+zb*q2+zc*q3,x), netv:listofvars(ha),lg:length(netv),
      if lg>0 then (ldf0(13),if ppowr(rb,rc,lg) then tp:70),
   if tp<=70 then go(ll2), vecdeg(za,zb,zc,x), tp:brchde2(rb,rc), 
   if tp=110 then (ldf0(2),svconf(rb,rc),return(ys)) else 
   if tp=120 then (ldf0(3),svhypr(rb,rc),return(ys)) else
   if tp=125 then (idvprep(2),go(ll1)),
   if tp=126 then (idvprep(1),go(ll1)),
   if tp=130 then return(ys),   if swsv=tdv or tp=990 then go(ll3),
  ll2,ldf0(5),neq:ivtrns(rb,rc,rr),if 10<tp and tp<14 then go(ll3),
   if neq#f  then (idvprep(1),go(ll1)),   if tp=0 then go(ll4),
  ll3, ldf0(6),neq:dvtrns(za,rb,rc,rr,ro,lam1),
   if neq#f  then (idvprep(2),go(ll1)),
  ll4, nosol(oeq),return(oeq))$

deqcoeff(exp):=block([],za:coeff(exp,dyx2,1),zb:coeff(exp,dyx1,1),
   zc:coeff(exp,y,1),rb:ratsimp(zb/za),rc:ratsimp(zc/za))$

brchde1(exp):=block([],
   if not freeof(sin,cos,tan,cot,exp)     then return(10),
   if not freeof(%e,exp)                  then return(20),
   if not freeof(log,exp)                 then return(30),
   if not freeof(sinh,cosh,tanh,coth,exp) then return(40),
   if not freeof(sn,cn,dn,pe,th,exp)      then return(50),
   if not freeof(sqrt(x),exp)             then return(60),
   return(990))$

brchde2(wb,wc):=block([ws,wi], 
   if nspt=1 then (if rk[1]=1 then return(110)),
   if nspt=2 then (if rk[1]<=0 and denom(rk[2])<3 and rk[2]<=1 then 
     (if spt[1]#0 then (axb:x-spt[1],return(91)) else return(110))),
   if nspt=3 and rk[1]=0 and rk[2]=0 then (if rk[3]=0 then
     (if nf=1 and maxacf=2 then (ldf0(14),if powr(wb,wc,2) then return(80)),
      if da=2 and db=1 and dc=0 then 
         (ldf0(15),tp:eqpt1(dum), if tp#f then return(tp)),
      if da=3 and db=2 and dc=1 then
         (ldf0(18),tp:eqpt12(wb,wc),if tp#f then return(tp))),
    if da=4 and db=3 and dc<=4 then
         (ldf0(17),tp:eqpt11(dum),if tp#f then return(tp)) else
    if rk[3]=0 then return(120)),
   ldf0(7),tp:brchde3(wb,wc),return(tp))$

tblspt(za,wb,wc,x,ksw):=block([k,n,wa,wf,wfi], k:0, n:hipow(za,x), 
   if n=0 then go(lf), pspt:1, maxacf:0, wf:facza(za), j:0,
   for i:1 thru nf do (pwf[i]:part(wf,i),if pwf[i]=x then j:i),
   if nf>1 and j>1 then (ws:pwf[1],pwf[1]:pwf[j],pwf[j]:ws),
   for i:1 thru nf do (wfi:solve(pwf[i],x),l1:length(wfi),
     lacf[i]:l1, l2:polord(za,pwf[i]), maxacf:max(maxacf,l2),
     for j:1 thru l1 do (k:k+1,spt[k]:part(wfi[j],2),
      if ksw#n then kindspt(wb,wc,k),
      if l1=1 then pspt:pspt*spt[k] else pspt:0)),
 lf,k:k+1,nspt:k,spt[k]:inf,if ksw#n then (kindspt(wb,wc,k),ro:rk[k]))$

facza(za):=block([w1,w2,w3,w4], w1:za,w2:diff(w1,x),w3:gcd(w1,w2), 
   w4:ratsimp(w1/w3), wf:factor(w4), if wf=w4 then nf:1 else nf:length(wf), 
   if nf=1 then wf:[wf], return(wf))$

vecdeg(za,zb,zc,x):=block([],da:hipow(za,x),db:hipow(zb,x),dc:hipow(zc,x))$

nosol(oeq):=block([], print("we cannot solve"), print(oeq))$

idvprep(i):=block([],print(str2,vtr),print(str3,neq), 
   if i=1 then eq:subst(x,t,neq) else eq:subst(y,u,neq), vprep(dum))$

vprep(dum):=block([],eq:stddeq(eq),print(str1,eq),kt:kt+1,aeq[kt]:eq,
   avtr[kt]:vtr)$

stddeq(oeq):=block([wexp,wr,deq],wexp:part(oeq,1),
   wr:part(oeq,2),deqcoeff(wexp),rr:ratsimp(wr/za), 
   za:calca(rb,rc,x), zb:ratsimp(rb*za), zc:ratsimp(rc*za), 
   deq:dyx2+rb*dyx1+rc*y=rr, return(deq))$

calca(wb,wc,x):=block([],w1:denom(wb),w2:denom(wc),za:lcm(w1,w2,x))$

lcm(p,q,x):=block([r,s],r:gcd(p,q,x),s:ratsimp(p*q/r),return(s))$

polord(exp,sexp):=block([we,n],we:exp,n:0,ll,we:ratsimp(we/sexp), 
   if denom(we)#1 then return(n), n:n+1,go(ll))$

kindspt(wb,wc,k):=block([bt,ct,b1,c1,dbt,dct],
   if spt[k]=inf then alpha:0 else alpha:spt[k],
   bt:ratsubst(1/t+alpha,x,wb), b1:ratsimp((2*t-bt)/t^2), 
   ct:ratsubst(1/t+alpha,x,wc), c1:ratsimp(ct/t^4),if spt[k]=inf then 
    (bt:subst(x,t,b1),ct:subst(x,t,c1),dbt:ratdeg(wb,x),dct:ratdeg(wc,x))
 else (bt:ratsimp(wb),ct:ratsimp(wc),  dbt:ratdeg(b1,t),dct:ratdeg(c1,t)),
   rk[k]:max(dbt,dct/2)+1, if rk[k]=0 and freeof(%i,alpha) 
   then regsingpt(x-alpha,bt,ct,alpha))$

ratdeg(rf,x):=block([wrf,nrf,dgnrf,drf,dgdrf,dgrf],if rf=0 then 
   return(-99), wrf:ratsimp(rf), nrf:num(wrf), dgnrf:hipow(nrf,x),
   drf:denom(wrf),dgdrf:hipow(drf,x), dgrf:ratsimp(dgnrf-dgdrf),
   return(dgrf))$

regsingpt(xa,wb,wc,alpha):=block([f1,g1,flx], remvalue(l),
   f1:ratsimp(xa*wb), g1:ratsimp(xa^2*wc), flx:l^2+(f1-1)*l+g1,
   cheq[k]:subst(alpha,x,flx),cflx[k]:flx)$

reqans1(dum):=block([],ans:read("continue?",str4))$

ldf0(i):=block([],if ldsw[0,i]#y then (ldsw[0,i]:y,
   if i= 2 then load(pconfl) else
   if i= 3 then load(phypgm) else
   if i= 5 then load(pivtr) else
   if i= 6 then load(pdvtr) else
   if i= 7 then load(pbrnch) else
   if i= 8 then load(pnonh) else
   if i= 9 then load(psode1) else
   if i=10 then load(pskelp) else
   if i=11 then load(pgptm) else
   if i=12 then load(premap) else
   if i=13 then load(ppow) else
   if i=14 then load(pnpow) else
   if i=15 then load(peqp1) else
   if i=16 then load(pstdsp) else
   if i=17 then load(peqp11) else
   if i=18 then load(peqp12) ))$

LODESAVE(FILENAME,[LODE_FUNCTIONS]):=
  IF STATUS(FEATURE,ITS)=TRUE
  THEN APPLY('FASSAVE,CONS(FILENAME,LODE_FUNCTIONS))
  ELSE BLOCK([PACKAGEFILE:TRUE],
	     IF STATUS(FEATURE,UNIX)=TRUE
	     THEN APPLY('SAVE,
			CONS(CONCAT(FIRST(FILENAME),".l"),lode_functions))
	     ELSE ERROR("unknown System"))$

lodesave([pmain,fasl],lode2,tobegin,nhstddeq,hlode2,brchde1,vecdeg,
   brchde2,deqcoeff,tblspt,facza,nosol,idvprep,vprep,stddeq,calca,
   lcm,polord,kindspt,ratdeg,regsingpt,reqans1,ldf0);

lodesave([LODSAV,FASL],LODESAVE)$