module odeprob! Brusselator with diffusion 200x200 grid
! This file bruss2h.f90 is part of the EPPEER package by B.A.Schmitt and R.Weiner
! available at www.mathematik.uni-marburg.de\~schmitt\peer
! It contains the Brusselator example from the documentation
!------------------------
  integer :: nprob = 80000, nparm = 1
  character(4) :: odename = "br2h"
  logical :: exsol = .false. ! solution unknown
contains

      pure subroutine fcn(t,u,fw,par)
!     Brusselator
      implicit none
      real(8), intent(in) :: t
      real(8), dimension(:), intent(in) :: u   ! exactly like rkstep
      real(8), dimension(:), intent(in) :: par
      real(8), dimension(:), intent(out) :: fw ! exactly like rkstep
      integer m, i,j,k,l,m2
      real*8 r,alf
      parameter (m=200, alf=1d-5)
      r=alf*(m-1)*(m-1)
      m2=m*m
      do 100 i=1,m2
         j=i+m2
         fw(i)=-4d0*u(i)
         fw(j)=-4d0*u(j)
 100  continue
      do 101 i=1,m2-m+1,m
         j=i+m2
         fw(i)=fw(i)+2d0*u(i+1)
         fw(j)=fw(j)+2d0*u(j+1)
         if (i.gt.1) then
            fw(i)=fw(i)+u(i-m)
            fw(j)=fw(j)+u(j-m)
         endif
         if (i.lt.m2-m) then
            fw(i)=fw(i)+u(i+m)
            fw(j)=fw(j)+u(j+m)
         endif
 101  continue
      do 102 i=1,m
         j=i+m2
         fw(i)=fw(i)+2d0*u(i+m)
         fw(j)=fw(j)+2d0*u(j+m)
         if (i.gt.1) then
            fw(i)=fw(i)+u(i-1)
            fw(j)=fw(j)+u(j-1)
         endif
         if (i.lt.m) then
            fw(i)=fw(i)+u(i+1)
            fw(j)=fw(j)+u(j+1)
         endif
 102  continue
      do 103 i=m,m2,m
         j=i+m2
         fw(i)=fw(i)+2d0*u(i-1)
         fw(j)=fw(j)+2d0*u(j-1)
         if (i.gt.m) then
            fw(i)=fw(i)+u(i-m)
            fw(j)=fw(j)+u(j-m)
         endif
         if (i.lt.m2) then
            fw(i)=fw(i)+u(i+m)
            fw(j)=fw(j)+u(j+m)
         endif
 103  continue
      do 104 i=m2-m+1,m2
         j=i+m2
         fw(i)=fw(i)+2d0*u(i-m)
         fw(j)=fw(j)+2d0*u(j-m)
         if (i.gt.m2-m+1) then
            fw(i)=fw(i)+u(i-1)
            fw(j)=fw(j)+u(j-1)
         endif
         if (i.lt.m2) then
            fw(i)=fw(i)+u(i+1)
            fw(j)=fw(j)+u(j+1)
         endif
 104  continue
      do 105 i=2,m-1
         do 105 k=1,m-2
          l=i+k*m
          j=l+m2
        fw(l)=fw(l)+u(l-1)+u(l+1)+u(l-m)+u(l+m)
        fw(j)=fw(j)+u(j-1)+u(j+1)+u(j+m)+u(j-m)
 105  continue
      do 106 i=1,m2
         j=i+m2
         fw(i)=r*fw(i)+1d0+u(i)*u(i)*u(j)-4d0*u(i)
         fw(j)=r*fw(j)+3d0*u(i)-u(i)*u(i)*u(j)
 106  continue
      return
      end subroutine

      subroutine inivals(t0,te,u,par)
      implicit none
!      Anfangswert Dgl
      integer n,m, i,j,k,l,m2
      real(8), dimension(:), intent(in) :: par
      real(8), dimension(:), intent(out) :: u
      real(8) t0,te,alf,dx,x,y
      parameter (m=200,alf=1d-5)
       n = 2*m*m
!      print "(a,i4)", "Brusselator n=", n
       m2=m*m
      t0 = 0D0
      te = 1D0
      dx = 1d0/(m-1)
      do 100 i=1,m
       do 100 k=1,m
       l=i+(k-1)*m
       y=(k-1)*dx
       u(l)=0.5d0+y
       j=l+m2
       x=(i-1)*dx
      u(j) = 1d0+5d0*x
 100  continue
      return
      end subroutine

!**********************************************************************

    subroutine solution(t,u,par)
      implicit none
	  integer n,i
	  real(8) t
      real(8), dimension(:), intent(out) :: u
      real(8), dimension(:), intent(in) :: par
!  unknown, approximate with exsol=.false.
      return
      end subroutine
end module odeprob