c-----------------------------------------------------------------------
c     these are the subroutines for the
c
c        Plate
c        ODE of dimension 80
c
c-----------------------------------------------------------------------
      subroutine prob(problm,neqn,ndisc,t,ijac,mljac,mujac)
      IMPLICIT double precision (A-H,O-Z)
      character*(*) problm
      integer neqn,ijac,mljac,mujac,ndisc
      double precision t(0:*)
      PARAMETER(MX=8,MY=5,MACHS1=2,MACHS2=4)
      COMMON /TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     +    OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT

      problm = 'plate'
      neqn   = 80
      ndisc  = 0
      t(0)   = 0d0
      t(1)   = 7.0d0
      ijac   = 1
      mljac  = neqn
      mujac  = neqn

C --- SET DEFAULT VALUES
      NX=MX
      NY=MY
      NACHS1=MACHS1
      NACHS2=MACHS2
      NXM1=NX-1
      NYm1=NY-1
      NDEMI=NX*NY
      OMEGA=1000.0D+0
      STIFFN=100.0D+0
      WEIGHT=200.0D+0
      DENOM=NX+1
      DELX=2.0D+0/DENOM
      USH4=1.0D+0/(DELX**4)
      FAC=STIFFN*USH4
      return
      end
c-----------------------------------------------------------------------
      subroutine init(neqn,y)
      integer neqn
      double precision y(neqn)

      integer i

      do 10 i=1,neqn
         y(i) = 0d0
   10 continue

      return
      end
c-----------------------------------------------------------------------
      SUBROUTINE feval(N, X, Y, DF, IERR, RPAR, IPAR)
      IMPLICIT real*8 (A-H,O-Z)
      DIMENSION Y(N), DF(N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
c ################ print
c      write(6,*)'  dans fplate, n,x=',n,x
c##############
C -------- LA BOUCLE -------
      DO 1 I=1,NX
      DO 1 J=1,NY
      K=I+NX*(J-1)
C -------- DERIVEE DEUXIEME ----
      DF(K)=Y(K+NDEMI)
C ------ POINT CENTRAL ---
      UC=16.D0*Y(K)
      IF(I.GT.1)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K-1)
      END IF
      IF(I.LT.NX)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K+1)
      END IF
      IF(J.GT.1)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K-NX)
      END IF
      IF(J.LT.NY)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K+NX)
      END IF
      IF(I.GT.1 .AND.J.GT.1 )UC=UC+2.D0*Y(K-NX-1)
      IF(I.LT.NX.AND.J.GT.1 )UC=UC+2.D0*Y(K-NX+1)
      IF(I.GT.1 .AND.J.LT.NY)UC=UC+2.D0*Y(K+NX-1)
      IF(I.LT.NX.AND.J.LT.NY)UC=UC+2.D0*Y(K+NX+1)
      IF(I.GT.2)UC=UC+Y(K-2)
      IF(I.LT.NXM1)UC=UC+Y(K+2)
      IF(J.GT.2)UC=UC+Y(K-2*NX)
      IF(J.LT.NYM1)UC=UC+Y(K+2*NX)
      IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN
         XI=I*DELX
         FORCE=EXP(-5.D0*(X-XI-2.D0)**2)+EXP(-5.D0*(X-XI-5.D0)**2)
      ELSE
         FORCE=0.D0
      END IF
      DF(K+NDEMI)=-OMEGA*Y(K+NDEMI)-FAC*UC+FORCE*WEIGHT
  1   CONTINUE
      RETURN
      END


c-----------------------------------------------------------------------
      subroutine jeval(N,X,Y,DFY,LDIM,IERR,IPAR,RPAR)
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),DFY(LDIM,N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
c ################ print
c      write(6,*)'  dans jac , n,x=',n,x
c##############
C------ METTRE A ZERO -------
      DO 1 I=1,N
      DO 1 J=1,N
  1   DFY(I,J)=0.D0
C -------- LA BOUCLE -------
      DO 2 I=1,NX
      DO 2 J=1,NY
      K=I+NX*(J-1)
C -------- DERIVEE DEUXIEME ----
      DFY(K,K+NDEMI)=1.D0
C ------ POINT CENTRAL ---
      DFY(K+NDEMI,K)=-FAC*16.D0
      IF(I.GT.1)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K-1)=FAC*8.D0
      END IF
      IF(I.LT.NX)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K+1)=FAC*8.D0
      END IF
      IF(J.GT.1)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K-NX)=FAC*8.D0
      END IF
      IF(J.LT.NY)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K+NX)=FAC*8.D0
      END IF
      IF(I.GT.1 .AND.J.GT.1 )DFY(K+NDEMI,K-NX-1)=-FAC*2.D0
      IF(I.LT.NX.AND.J.GT.1 )DFY(K+NDEMI,K-NX+1)=-FAC*2.D0
      IF(I.GT.1 .AND.J.LT.NY)DFY(K+NDEMI,K+NX-1)=-FAC*2.D0
      IF(I.LT.NX.AND.J.LT.NY)DFY(K+NDEMI,K+NX+1)=-FAC*2.D0
      IF(I.GT.2)DFY(K+NDEMI,K-2)=-FAC
      IF(I.LT.NXM1)DFY(K+NDEMI,K+2)=-FAC
      IF(J.GT.2)DFY(K+NDEMI,K-2*NX)=-FAC
      IF(J.LT.NYM1)DFY(K+NDEMI,K+2*NX)=-FAC
      DFY(K+NDEMI,K+NDEMI)= -OMEGA
  2   CONTINUE
      RETURN
      END


c-----------------------------------------------------------------------
      subroutine solut(neqn,true)
      integer neqn
      double precision true(neqn)
c
c

         true(1)=0.490143813851336d-3
         true(2)=0.980081485560611d-3
         true(3)=0.1462893811482190d-2
         true(4)=0.1915822464411935d-2
         true(5)=0.2285152533727002d-2
         true(6)=0.2461353376688549d-2
         true(7)=0.2254597413097122d-2
         true(8)=0.1438312591933600d-2
         true(9)=0.849025149228402d-3
         true(10)=0.1697885005625757d-2
         true(11)=0.2535239886068847d-2
         true(12)=0.3323989552181772d-2
         true(13)=0.3977902193560667d-2
         true(14)=0.4320231736082990d-2
         true(15)=0.4025679955083897d-2
         true(16)=0.2643206356123840d-2
         true(17)=0.980287627702671d-3
         true(18)=0.1960162971121222d-2
         true(19)=0.2925787622964379d-2
         true(20)=0.3831644928823870d-2
         true(21)=0.4570305067454005d-2
         true(22)=0.4922706753377098d-2
         true(23)=0.4509194826194244d-2
         true(24)=0.2876625183867201d-2
         true(25)=0.849025149228402d-3
         true(26)= 0.1697885005625757d-2
         true(27)=0.2535239886068847d-2
         true(28)=0.3323989552181772d-2
         true(29)=0.3977902193560667d-2
         true(30)=0.4320231736082990d-2
         true(31)=0.4025679955083897d-2
         true(32)=0.2643206356123840d-2
         true(33)=0.490143813851336d-3
         true(34)=0.980081485560611d-3
         true(35)=0.1462893811482190d-2
         true(36)=0.1915822464411935d-2
         true(37)=0.2285152533727002d-2
         true(38)=0.2461353376688549d-2
         true(39)=0.2254597413097122d-2
         true(40)=0.1438312591933600d-2
         true(41)=-0.1177590304545409d-2
         true(42)=-0.2409005827992214d-2
         true(43)=-0.3722140831656533d-2
         true(44)=-0.5078780056048207d-2
         true(45)=-0.6302661811097914d-2
         true(46)=-0.6973399942926759d-2
         true(47)=-0.6394575120415784d-2
         true(48)=-0.3960464551310118d-2
         true(49)=-0.2040148244040460d-2
         true(50)=-0.4174829877953482d-2
         true(51)=-0.6456510337516159d-2
         true(52)=-0.8832503276738242d-2
         true(53)=-0.11029624807177369d-1
         true(54)=-0.12352389570141255d-1
         true(55)=-0.11524177328690540d-1
         true(56)=-0.7253301886026949d-2
         true(57)=-0.2355180609090818d-2
         true(58)=-0.4818011655984428d-2
         true(59)=-0.7444281663313065d-2
         true(60)=-0.10157560112096413d-1
         true(61)=-0.12605323622195828d-1
         true(62)=-0.13946799885853517d-1
         true(63)=-0.12789150240831569d-1
         true(64)=-0.7920929102620235d-2
         true(65)=-0.2040148244040460d-2
         true(66)=-0.4174829877953482d-2
         true(67)=-0.6456510337516159d-2
         true(68)=-0.8832503276738242d-2
         true(69)=-0.11029624807177369d-1
         true(70)=-0.12352389570141255d-1
         true(71)=-0.11524177328690540d-1
         true(72)=-0.7253301886026949d-2
         true(73)=-0.1177590304545409d-2
         true(74)=-0.2409005827992214d-2
         true(75)=-0.3722140831656533d-2
         true(76)=-0.5078780056048207d-2
         true(77)=-0.6302661811097914d-2
         true(78)=-0.6973399942926759d-2
         true(79)=-0.6394575120415784d-2
         true(80)=-0.3960464551310118d-2

      return
      end