/*  THIS LITTLE PACKAGE SOLVES FIRST ORDER LINEAR
     ORDINARY DIFFERENTIAL EQUATIONS SUBJECT TO A 
     BOUNDARY CONDITION (B.C.) AT AN INITIAL POINT
       THE CALLING PROCEDURE IS
          IVPSOL(DIFFEQ,Y,X,A,BCEQ);
     WHERE DIFFEQ IS THE DIFFERENTIAL EQUATION TO BE SOLVED
            (E.G.:  X*'DIFF(Y,X)+2*Y=1)
           Y IS THE DEPENDENT VARIABLE
           X IS THE INDEPENDENT VARIABLE
           A IS THE POINT AT WHICH THE B.C. IS APPLIED
           BCEQ IS THE B.C. EQUATION (E.G.: Y=2)
*/
SOLDIFF(EQ1,Y,X,A):=
     BLOCK([EQ2,A1,A2,A3,A4,A5],
        EQ2: LHS(EQ1)-RHS(EQ1),
        A1:  RATCOEF(EQ2,'DIFF(Y,X)),
        A2:  RATCOEF(EQ2,Y),
        A3:  EQ2-A1*'DIFF(Y,X)-A2*Y,
        A4:  INTEGRATE(A2/A1,X),
        A5:  INTEGRATE(%E**A4*A3/A1,X),
           'CONST1/%E**A4-A5/%E**A4)$
IVPSOL(DIFFEQ,Y,X,A,BCEQ):=
     BLOCK([Z,DERY,DYA,YA,CO1],
        Z:   SOLDIFF(DIFFEQ,Y,X,A),
        DERY:DIFF(Z,X,1),
        DYA: SUBST(A,X,DERY),
        YA:  SUBST(A,X,Z),
        CO1: SOLVE(SUBST(YA,Y,SUBST(DYA,'DIFF(Y,X),BCEQ)),'CONST1),
        IF length(co1)>1 /* LISTP(CO1) */ THEN
           [PRINT("SOLUTIONS ARE NOT UNIQUE, POSSIBLE SOLUTIONS ARE:
             "),
            MAP(LAMBDA([L1],apply('DISPLAY,[L1])),CO1), RETURN(Z)] 
           ELSE   if co1=[] then return(print("No solution for given initial
condition"),z) else SUBST(RHS(CO1[1]),'CONST1,Z))$
/* FIND RIGHT EIGENVECTORS WITH OTHER COMPONENT = 1
          I.E. X IS SOLUTION OF M*X=MU*X    */
EVEC(M,MU,MODES):=block([equations,unknowns,x],
(FOR L:1 THRU MODES DO[
     EQUATIONS:[],
     UNKNOWNS:[],
/*     KILL(X), */
     X[L]:1,
     FOR I:1 THRU MODES DO IF NOT (I=L) THEN [
         UNKNOWNS:CONS(X[I],UNKNOWNS),
         EQUATIONS:CONS(SUM(M[I,J]*X[J],J,1,MODES)
              =MU[L]*X[I],EQUATIONS)],
     PRINT(SOLVE(EQUATIONS,UNKNOWNS))]))$
/* CONSTRUCT 3 X 3 MATRIX WITH DESIRED EIGENVALUES */

CONSMAT3(L1,L2,L3):=BLOCK([A,B,MAT3,CP,CPM,SOL1,L],
MAT3:MATRIX([1,2,3],[B,3,A],[1,1,L1+L2+L3-4]),
CP:EXPAND(-(L-L1)*(L-L2)*(L-L3)),
CPM:CHARPOLY(MAT3,L),GLOBALSOLVE:TRUE,
SOL1:SOLVE([RATCOEF(CP,L)=RATCOEF(CPM,L),
RATCOEF(CP,L,0)=RATCOEF(CPM,L,0)],[A,B]),
RETURN(apply('EV,[MAT3])))$