define_variable(g,'g,any)$

/*THIS BLOCK CHECKS FOR A POLYNOMIAL IN N*/
POLYP(G,N):=BLOCK([D,F,C],
     G:RATEXPAND(G), IF FREEOF(N,G) THEN RETURN(TRUE),
     D:HIPOW(G,N), F:TRUE,
     FOR I:D STEP -1 THRU 0 DO
     (C:COEFF(G,N,I), IF NOT(FREEOF(N,C)) THEN F:FALSE,
      G:RATEXPAND(G-C*N^I)),
     RETURN(IS(G=0 AND F)))$
   


/*THIS BLOCK CHECKS FOR A CONSTANT TO A POLYNOMIAL POWER*/
POLYINN(X,N):=BLOCK([B,E],
     IF INPART(X,0)="*" THEN
     RETURN(POLYINN(INPART(G,1),N) AND POLYINN(INPART(G,2),N)),
     IF INPART(X,0)#"^" THEN RETURN(FALSE),
     B:INPART(X,1),
     E:INPART(X,2),
     IF NOT FREEOF(N,B) THEN RETURN(FALSE),
     RETURN(POLYP(E,N)))$



/*THIS BLOCK IMPLEMENTS THE CHARACTERISTIC EQUATION METHOD*/
CHAR(E,G,U,N,K,IV):=BLOCK([GENSOL,HOMSOL,PARSOL,LOS,MULTIPLICITIES,
H,V,L,SS,DISPFLAG],
LOCAL(A,AA,B,R,M),
DISPFLAG:FALSE,
  FOR I:0 THRU K DO
    AA[I]:COEFF(E,U(N+K-I)),
    H:0,
  FOR I:0 THRU K DO
    H:H+AA[I]*U(N+K-I),
  IF H#E THEN RETURN("ERRONEOUS INPUT"),

  FOR I:0 THRU K DO
    H:SUBST(U^(K-I),U(N+K-I),H),

  MULTIPLICITIES:TRUE,
  LOS:SOLVE(H,U),
  FOR I:1 THRU LENGTH(LOS) DO
    (R[I]:LOS[I], R[I]:RHS(EV(R[I])),
     M[I]:MULTIPLICITIES[I]),
  HOMSOL:
    SUM(SUM(A[I,J]*N^(M[I]-J),J,1,M[I])*R[I]^N,I,1,LENGTH(LOS)),


  IF G=0 THEN
     (V:[ ],
      FOR I:1 THRU LENGTH(LOS) DO
      FOR J:1 THRU M[I] DO V:CONS(A[I,J],V),
         L:[ ],
      FOR Q:0 THRU K-1 DO L:CONS(SUBST(Q,N,HOMSOL)=U(Q),L),
      SS:EV(SOLVE(L,V),IV),
      RETURN(U(N)=(EV(HOMSOL,SS))))

  ELSE IF POLYP(G,N) = TRUE THEN
     (G:RATEXPAND(G), PARSOL:SUM(B[J]*N^J,J,0,HIPOW(G,N)),
      FOR J:0 THRU K DO
       (L:0, V:E,
         FOR I:0 THRU K DO
         (L:RATEXPAND(SUBST(N+K-I,N,B[J]*N^J)),
          V:RATEXPAND(SUBST(L,U(N+K-I),V))),
          V:RATSIMP(V),
           IF V#0 THEN RETURN(V) ELSE PARSOL:N*PARSOL),
         V:E,
      FOR I:0 THRU K DO (L:RATEXPAND(SUBST(N+K-I,N,PARSOL)),
         V:RATEXPAND(SUBST(L,U(N+K-I),V))),
         L:[ ],
      FOR I:0 THRU HIPOW(PARSOL,N) DO
         L:CONS(COEFF(V=G,N,I),L),
         V:[ ],
      FOR J:0 THRU HIPOW(PARSOL,N) DO
         V:CONS(B[J],V),
      SS:SOLVE(L,V),
      PARSOL:EV(PARSOL,SS))

  ELSE IF POLYINN(G,N) = TRUE THEN
     (PARSOL:B1*G,
      FOR J:0 THRU K DO
       (L:0, V:E,
         FOR I:0 THRU K DO
         (L:SUBST(N+K-I,N,PARSOL), V:SUBST(L,U(N+K-I),V)),
          V:RATSIMP(V),
          IF V#0 THEN RETURN(V) ELSE PARSOL:N*PARSOL),
      SS:SOLVE(V=G,B1),
      PARSOL:EV(PARSOL,SS))

  ELSE IF INPART(G,0)=SIN OR INPART(G,0) = COS THEN
     (PARSOL:B[1]*SIN(INPART(G,1)) + B[2]*COS(INPART(G,1)),
      FOR J:0 THRU K DO
       (L:0, V:E,
         FOR I:0 THRU K DO
         (L:EXPAND(SUBST(N+K-I,N,PARSOL)),
          V:EXPAND(SUBST(L,U(N+K-I),V))),
      V:TRIGEXPAND(V),
      IF V#0 THEN RETURN(V) ELSE PARSOL:N*PARSOL),
         V:E,
      FOR I:0 THRU K DO(L:EXPAND(SUBST(N+K-I,N,PARSOL)),
         V:EXPAND(SUBST(L,U(N+K-I),V))),
      V:TRIGEXPAND(V),
         L:[ ],
      LT:[SIN(INPART(G,1)),COS(INPART(G,1))],
      FOR JJ:1 THRU 2 DO
         L:CONS(COEFF(V=G,LT[JJ]),L),
         V:[ ],
      FOR J:1 THRU 2 DO
         V:CONS(B[J],V),
      SS:SOLVE(L,V),
      PARSOL:EV(PARSOL,SS))


  ELSE RETURN("CAN'T BE SOLVED IN CLOSED FORM BY PROGRAM"),

  GENSOL:HOMSOL + PARSOL,
     V:[ ],
  FOR I:1 THRU LENGTH(LOS) DO
  FOR J:1 THRU M[I] DO V:CONS(A[I,J],V),
     L:[ ],
  FOR Q:0 THRU K-1 DO
  L:CONS(SUBST(Q,N,GENSOL)=U(Q),L),
  SS:EV(SOLVE(L,V),IV),
  RETURN(U(N)=(EV(GENSOL,SS))))$




/*THIS BLOCK IMPLEMENTS THE GENERATING FUNCTION METHOD*/
GENF(E,G,U,N,K,IV):=BLOCK([MULTIPLICITIES,L,V,SS,VV,LOS,
NR,F,SOL,P,DISPFLAG],
LOCAL(A,AA,B),
DISPFLAG:FALSE,
  FOR I:0 THRU K DO
    AA[I]:COEFF(E,U(N+K-I)),
    H:0,
  FOR I:0 THRU K DO
    H:H+AA[I]*U(N+K-I),
  IF H#E THEN RETURN("ERRONEOUS INPUT"),

     L:E,
  FOR I:0 THRU K DO
    L:SUBST((F-SUM(U(J)*X^J,J,0,K-I-1))*X^I,U(N+K-I),L),

  IF G=0 THEN
     (S:SOLVE(L,F),
      F:EV(F,S))

  ELSE IF POLYP(G,N) = TRUE THEN
     (G:RATEXPAND(G),
      V:SUBST(X^K/(1-X)*COEFF(G,N,0),COEFF(G,N,0),G),
      VV:RATSIMP(DIFF(1/(1-X),X)),
      FOR I:1 THRU HIPOW(G,N) DO
        (V:SUBST(X^K*X*VV*COEFF(G,N,I),COEFF(G,N,I)*N^I,V),
         VV:RATSIMP(DIFF(X*VV,X))),
      V:RATSIMP(V),
      SS:SOLVE(L=V,F),
      F:EV(F,SS))

  ELSE IF POLYINN(G,N) = TRUE AND HIPOW(INPART(G,2),N) < 2 THEN
     (G1:(X^K)*(INPART(G,1)^COEFF(INPART(G,2),N,0)),
      G2:1 - X*(INPART(G,1)^COEFF(INPART(G,2),N,1)),
      V:RATSIMP(G1/G2),
      SS:SOLVE(L=V,F),
      F:EV(F,SS))

  ELSE RETURN("CAN'T BE SOLVED IN CLOSED FORM BY PROGRAM"),


  MULTIPLICITIES:TRUE,
  LOS:SOLVE(NEWRAT(F),X),
  FOR I:1 THRU LENGTH(LOS) DO
     (R[I]:LOS[I], R[I]:RHS(EV(R[I])),
      M[I]:MULTIPLICITIES[I]),

     V:[ ],
  B:PRODUCT((1-R[I]*X)^M[I],I,1,LENGTH(LOS)),
  FOR I:1 THRU LENGTH(LOS) DO
  FOR J:1 THRU M[I] DO
     (P[I,J]:B*A[I,J]/((1-R[I]*X)^J),  V:CONS(A[I,J],V)),
  P:SUM(SUM(P[I,J],J,1,M[I]),I,1,LENGTH(LOS)),

     L:[ ],
  NF:RATEXPAND(NUM(F)/ABS(COEFF(DENOM(F),X,0))), P:RATEXPAND(P),
  FOR I:0 THRU HIPOW(RATEXPAND(B),X)-1  DO
     L:CONS(COEFF(NF=P,X,I),L),
  SSS:EV(SOLVE(L,V),IV),

  SOL:SUM(SUM(A[I,J]*COEFF(DENOM(F),X,0)/ABS(COEFF(DENOM(F),X,0))*
     BINOMIAL(J+N-1,N)*R[I]^N,J,1,M[I]),I,1,LENGTH(LOS)),

  RETURN(U(N)=(EV(SOL,SSS))))$


/*THIS BLOCK FINDS THE NEW POLYNOMIAL ASSOCIATED TO F*/
NEWRAT(F):=BLOCK([HD,CP,DP],
     HD:HIPOW(DENOM(F),X),
     CP:COEFF(DENOM(F),X,HD),
     DP:SUM((COEFF(DENOM(F),X,I))/CP*X^I,I,0,HD),
     RETURN(SUM(COEFF(DP,X,HD-I)*X^I,I,0,HD)))$





/*THIS BLOCK IMPLEMENTS THE VARIABLE COEFFICIENT METHOD*/
VARC1(E,G,U,N,K,IV):=BLOCK([V,VV,EQ,Y,CAUCHYSUM,FINSOL,SERSOL,DISPFLAG],
LOCAL(A,B),DISPFLAG:FALSE,
  FOR I:0 THRU K DO
    (A[I]:COEFF(E,U(N+I)),
     A[I]:RATEXPAND(A[I]),
     IF POLYP(A[I],N)=FALSE THEN RETURN("CAN'T DO IT")),
  IF K=2 AND (B:BESSELCHECK(E,K) # FALSE) THEN RETURN(B),
  V:RATEXPAND(E),
  FOR I:K STEP -1 THRU 0 DO
     FOR J:HIPOW(A[I],N) STEP -1 THRU 0 DO
       (V:RATSUBST(X^J*'DIFF(Y,X,I+J),N^J*U(N+I),V),
        V:RATEXPAND(V)),
  V:RATSUBST(Y,'DIFF(Y,X,0),V),
  V:RATEXPAND(V),
  IF POLYP(G,N) = TRUE THEN
    (G:RATEXPAND(G), VV:G,
     FOR I:0 THRU HIPOW(G,N) DO
        VV:SUBST(X^I,N^I,VV),
        VV:%E^X*VV)
  ELSE RETURN("CAN'T DO IT"),

     EQ:V-VV,
  DEPENDENCIES(Y(X)),
  SOL:ODE2(EQ=0,Y,X),
        IF K=1 THEN FINSOL:IC1(SOL,X=0,Y=EV(U(0),IV))
        ELSE IF K=2 THEN FINSOL:IC2(SOL,X=0,Y=EV(U(0),IV),'DIFF(Y,X)=EV(U(1),IV))
  ELSE RETURN("O.D.E. CAN'T BE SOLVED AT PRESENT BY MACSYMA"),
     CAUCHYSUM:TRUE,
  SERSOL:POWERSERIES(RHS(FINSOL),X,0), SERSOL:EXPAND(SERSOL),
  B:INPART(SERSOL,1),
  B:EV(B,X=1),
  IF ATOM(B)=FALSE THEN B:SUBSTPART(N,B,4),
  RETURN(U(N)=((N!)*B)))$


/*THIS BLOCK CHECKS FOR A BESSEL RECURRENCE RELATION*/

BESSELCHECK(E,K):=BLOCK([A,ANS],
LOCAL(A),
     FOR I:0 THRU K DO
     (A[I]:COEFF(E,U(N+I)),
      A[I]:RATEXPAND(A[I])),
     IF NOT(INTEGERP(A[0])) THEN RETURN(FALSE),
     IF NOT(INTEGERP(EV(A[1],N=0))) THEN RETURN(FALSE),
     IF NOT(HIPOW(A[1],N)=1) THEN RETURN(FALSE),
     IF NOT(INTEGERP(COEFF(A[1],N,1))) THEN RETURN(FALSE),
     IF NOT(A[2]=1) THEN RETURN(FALSE),
     ANS:"A LINEAR COMBINATION OF BESSEL FUNCTIONS",
     /*EXACT DETAILS ARE OF NO SIGNIFICANCE,SINCE WE ARE MERELY
       DEMONSTRATING THE FEASIBILITY OF THIS APPROACH*/
     RETURN(ANS))$





/*THIS BLOCK IMPLEMENTS THE FIRST ORDER METHOD*/
VARC2(E,G,U,N,K,IV):=BLOCK([H,P,V,C,SOL],
LOCAL(AP,P),
     P:(-1)*COEFF(E,U(N))/COEFF(E,U(N+1)),
     V:G/COEFF(E,U(N+1)),
     S[J]:SUBST(J,N,P),
     S[I]:SUBST(I,N,P),
     P[N]:PRODUCT(S[I],I,1,N-1),
     H[I]:SUBST(I,N,V)/PRODUCT(S[J],J,1,I-1),
     V1:SUM(H[I],I,0,N),
     AP:EV(U(0)-SUBST(0,N,V),IV),
     RETURN(U(N)=AP*P[N]+P[N]*V1))$