C THIS PROGRAM COMPUTES THE SOLUTION OF THE EXAMPLE PROBLEM GIVEN
C IN THE DOCUMENTATION FOR SUBROUTINE MUSL.
C
DOUBLE PRECISION A,B,MA(3,3),MB(3,3),BCV(3),AMP,ER(5),TI(15),
1 X(3,15),U(6,15),Q(3,3,15),D(3,15),PHIREC(6,15),W(42),
2 EXSOL,AE
INTEGER IW(9)
EXTERNAL FLIN,FDIF
C
C SETTING OF THE INPUT PARAMETERS
C
N = 3
IHOM = 1
ER(1) = 1.D-11
ER(2) = 1.D-6
ER(3) = 1.1D-16
NRTI = 10
NTI = 15
NU = 6
LW = 42
LIW = 9
A = 0.D0
B = 6.D0
AMP = 0.D0
C
C SETTING THE BC MATRICES MA AND MB
C
DO 1100 I = 1 , N
DO 1000 J = 1 , N
MA(I,J) = 0.D0
MB(I,J) = 0.D0
1000 CONTINUE
MA(I,I) = 1.D0
MB(I,I) = 1.D0
1100 CONTINUE
C
C SETTING THE BC VECTOR BCV
C
BCV(1) = 1.D0 + DEXP(6.D0)
BCV(2) = BCV(1)
BCV(3) = BCV(1)
C
C CALL MUSL
C
CALL MUSL(FLIN,FDIF,N,IHOM,A,B,MA,MB,BCV,AMP,ER,NRTI,TI,NTI,
1 X,U,NU,Q,D,KPART,PHIREC,W,LW,IW,LIW,IERROR)
IF ((IERROR.NE.0).AND.(IERROR.NE.200).AND.(IERROR.NE.213).AND.
1 (IERROR.NE.240)) GOTO 5000
C
C COMPUTATION OF THE ABSOLUTE ERROR IN THE SOLUTION AND WRITING
C OF THE SOLUTION AT THE OUTPUTPOINTS
C
WRITE(*,200)
WRITE(*,190) ER(4),ER(5)
WRITE(*,210)
WRITE(*,200)
DO 1500 K = 1 , NRTI
EXSOL = DEXP(TI(K))
AE = EXSOL - X(1,K)
WRITE(*,220) K,TI(K),X(1,K),EXSOL,AE
DO 1300 I = 2 , N
AE = EXSOL - X(I,K)
WRITE(*,230) X(I,K),EXSOL,AE
1300 CONTINUE
1500 CONTINUE
STOP
5000 WRITE(*,300) IERROR
STOP
C
190 FORMAT(' CONDITION NUMBER = ',D10.3,/,
1 ' AMPLIFICATION FACTOR = ',D10.3,/)
200 FORMAT(' ')
210 FORMAT(' I ',6X,'T',8X,'APPROX. SOL.',9X,'EXACT SOL.',8X,
1 'ABS. ERROR')
220 FORMAT(' ',I3,3X,F7.4,3(3X,D16.9))
230 FORMAT(' ',13X,3(3X,D16.9))
300 FORMAT(' TERMINAL ERROR IN MUSL: IERROR = ',I4)
C
END
C
SUBROUTINE FLIN(T,Y,F)
C ----------------------
C
DOUBLE PRECISION T,Y(3),F(3)
DOUBLE PRECISION TI,SI,CO
C
TI = 2.D0 * T
SI = 2.D0 * DSIN(TI)
CO = 2.D0 * DCOS(TI)
F(1) = (1.D0 - CO) * Y(1) + (1.D0 + SI) * Y(3)
F(2) = 2.D0 * Y(2)
F(3) = (-1.D0 + SI) * Y(1) + (1.D0 + CO) * Y(3)
C
RETURN
C END OF FLIN
END
C
SUBROUTINE FDIF(T,Y,F)
C ----------------------
C
DOUBLE PRECISION T,Y(3),F(3)
DOUBLE PRECISION TI,SI,CO
C
CALL FLIN(T,Y,F)
TI = 2.D0 * T
SI = 2.D0 * DSIN(TI)
CO = 2.D0 * DCOS(TI)
TI = DEXP(T)
F(1) = F(1) + (-1.D0 + CO - SI)*TI
F(2) = F(2) - TI
F(3) = F(3) + (1.D0 - CO - SI)*TI
C
RETURN
C END OF FDIF
END
.