c       This is the problem definition for Brain Tumor Model.
c       PDE: D_x*u_x + D*u_xx, u_x(a,t)=0, u_x(b,t)=0, 
c            u(x,0) = (exp(-((x - xi) * (x - xi))/(eta*eta)))/(eta*sqrt(pi))
c       This code uses loops to produce multiple independent copies of
c       the problem in order to artificially increase the computation
c       required, if so desired. Loop bounded by npde.
c-----------------------------------------------------------------------
      subroutine f(t, x, u, ux, uxx, fval, npde)
c-----------------------------------------------------------------------
c purpose:
c       this subroutine defines the right hand side vector of the
c       npde dimensional parabolic partial differential equation
c                        ut = f(t, x, u, ux, uxx).
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 npde
c                               the number of pdes in the system.
c
        double precision        t
c                               the current time coordinate.
c
        double precision        x
c                               the current spatial coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
        double precision        uxx(npde)
c                               uxx(1:npde) is the approximation of the
c                               second spatial derivative of the
c                               solution at the point (t,x).
c
c output:
        double precision        fval(npde)
c                               fval(1:npde) is the right hand side
c                               vector f(t, x, u, ux, uxx) of the pde.
c-----------------------------------------------------------------------
        double precision a, b, w1, w2, gamma, xi, eta
	common /tumor/ a, b, w1, w2, gamma, xi, eta
c-----------------------------------------------------------------------
c
c local variables
        double precision        l, D, Dx, temp
c-----------------------------------------------------------------------
c loop indices:
        integer                 i
c-----------------------------------------------------------------------
c
c     assign fval(1:npde) according to the right hand side of the pde
c     in terms of u(1:npde), ux(1:npde), uxx(1:npde).
c
	  l = 30.0D0
	  D = (1.0D0/(exp(l * (w1 - x)) + 1.0D0)) - 1.0D0
	  D = (D + (1.0D0/(exp(l * (x - w2)) + 1.0D0))) * (1.0D0 - gamma)
	  D = D + gamma
	  Dx = (exp(l *(w1 - x))/((exp(l * (w1 - x)) + 1.0D0) ** 2.0D0))
	  temp = (exp(l *(x - w2))/((exp(l * (x - w2)) + 1.0D0) ** 2.0D0))
	  Dx = (Dx - temp) * (1.0D0 - gamma) * l

      if (isnan(Dx)) then
        Dx = 0 
      endif

      do i = 1, npde
         fval(i) = D*uxx(i) + Dx*ux(i)
      end do
c
      return
      end
c-----------------------------------------------------------------------
      subroutine derivf(t, x, u, ux, uxx, dfdu, dfdux, dfduxx, npde)
c-----------------------------------------------------------------------
c purpose:
c       this subroutine is used to define the information about the
c       pde required to form the analytic jacobian matrix for the dae
c       or ode system. assuming the pde is of the form
c                        ut = f(t, x, u, ux, uxx)
c       this routine returns the jacobians d(f)/d(u), d(f)/d(ux), and
c       d(f)/d(uxx).
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 npde
c                               the number of pdes in the system.
c
        double precision        t
c                               the current time coordinate.
c
        double precision        x
c                               the current spatial coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
        double precision        uxx(npde)
c                               uxx(1:npde) is the approximation of the
c                               second spatial derivative of the
c                               solution at the point (t,x).
c
c output:
        double precision        dfdu(npde,npde)
c                               dfdu(i,j) is the partial derivative
c                               of the i-th component of the vector f
c                               with respect to the j-th component
c                               of the unknown function u.
c
        double precision        dfdux(npde,npde)
c                               dfdux(i,j) is the partial derivative
c                               of the i-th component of the vector f
c                               with respect to the j-th component
c                               of the spatial derivative of the
c                               unknown function u.
c
        double precision        dfduxx(npde,npde)
c                               dfduxx(i,j) is the partial derivative
c                               of the i-th component of the vector f
c                               with respect to the j-th component
c                               of the second spatial derivative of the
c                               unknown function u.
c-----------------------------------------------------------------------
      double precision a, b, w1, w2, gamma, xi, eta
	  common /tumor/ a, b, w1, w2, gamma, xi, eta
c-----------------------------------------------------------------------
c
c local variables
        double precision        l, D, Dx, temp
c-----------------------------------------------------------------------
c loop indices:
        integer                 i, j
c-----------------------------------------------------------------------
c
c     assign dfdu(1:npde,1:npde), dfdux(1:npde,1:npde), and
c     dfduxx(1:npde,1:npde) according to the right hand side of the pde
c     in terms of u(1:npde), ux(1:npde), uxx(1:npde).
c
	  l = 30.0D0
	  D = (1.0D0/(exp(l * (w1 - x)) + 1.0D0)) - 1.0D0
	  D = (D + (1.0D0/(exp(l * (x - w2)) + 1.0D0))) * (1.0D0 - gamma)
	  D = D + gamma
	  Dx = (exp(l *(w1 - x))/((exp(l * (w1 - x)) + 1.0D0) ** 2.0D0))
	  temp = (exp(l *(x - w2))/((exp(l * (x - w2)) + 1.0D0) ** 2.0D0))
	  Dx = (Dx - temp) * (1.0D0 - gamma) * l

      if (isnan(Dx)) then
        Dx = 0 
      endif
	  
      do i = 1, npde
         do j = 1, npde
            dfdu(i,j) = 0.d0
            dfdux(i,j) = 0.d0
            dfduxx(i,j) = 0.d0
         end do
		 dfdux(i,i) = Dx
         dfduxx(i,i) = D
      end do
c
      return
      end
c-----------------------------------------------------------------------
      subroutine bndxa(t, u, ux, bval, npde)
c-----------------------------------------------------------------------
c purpose:
c       the subroutine is used to define the boundary conditions at the
c       left spatial end point x = xa.
c                           b(t, u, ux) = 0
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 npde
c                               the number of pdes in the system.
c
        double precision        t
c                               the current time coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,xa).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,xa).
c
c output:
        double precision        bval(npde)
c                               bval(1:npde) is the boundary condition
c                               at the left boundary point.
c-----------------------------------------------------------------------
c loop indices:
        integer                 i
c-----------------------------------------------------------------------
      do i = 1, npde
         bval(i) = ux(i)
      end do
c
      return
      end
c-----------------------------------------------------------------------
      subroutine bndxb(t, u, ux, bval, npde)
c-----------------------------------------------------------------------
c purpose:
c       the subroutine is used to define the boundary conditions at the
c       right spatial end point x = xb.
c                           b(t, u, ux) = 0
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 npde
c                               the number of pdes in the system.
c
        double precision        t
c                               the current time coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,xb).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,xb).
c
c output:
        double precision        bval(npde)
c                               bval(1:npde) is the boundary condition
c                               at the right boundary point.
c-----------------------------------------------------------------------
c loop indices:
        integer                 i
c-----------------------------------------------------------------------
      do i = 1, npde
         bval(i) = ux(i)
      end do
c
      return
      end
c-----------------------------------------------------------------------
      subroutine difbxa(t, u, ux, dbdu, dbdux, dbdt, npde)
c-----------------------------------------------------------------------
c purpose:
c       the subroutine is used to define the differentiated boundary
c       conditions at the left spatial end point x = xa. for the
c       boundary condition equation
c                              b(t, u, ux) = 0
c       the partial derivatives db/du, db/dux, and db/dt are supplied
c       by this routine.
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 npde
c                               the number of pdes in the system.
c
        double precision        t
c                               the current time coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
c output:
        double precision        dbdu(npde,npde)
c                               dbdu(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the unknown function u.
c
        double precision        dbdux(npde,npde)
c                               dbdux(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the spatial derivative of the
c                               unknown function u.
c
        double precision        dbdt(npde)
c                               dbdt(i) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to time t.
c-----------------------------------------------------------------------
c loop indices:
        integer                 i, j
c-----------------------------------------------------------------------
c
c     assign dbdu(1:npde,1:npde), dbdu(1:npde,1:npde), and dbdt(1:npde)
c     according to the right boundary condition equation in terms of
c     u(1:npde), ux(1:npde), uxx(1:npde).
c
      do i = 1, npde
         do j = 1, npde
            dbdu(i,j) = 0.0d0
            dbdux(i,j) = 0.0d0
         end do
         dbdt(i) = 0.d0
		 dbdux(i,i) = 1.0d0
      end do
c
      return
      end
c-----------------------------------------------------------------------
      subroutine difbxb(t, u, ux, dbdu, dbdux, dbdt, npde)
c-----------------------------------------------------------------------
c purpose:
c       the subroutine is used to define the differentiated boundary
c       conditions at the right spatial end point 1 = xb. for the
c       boundary condition equation
c                              b(t, u, ux) = 0
c       the partial derivatives db/du, db/dux, and db/dt are supplied
c       by this routine.
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 npde
c                               the number of pdes in the system.
c
        double precision        t
c                               the current time coordinate.
c
        double precision        u(npde)
c                               u(1:npde) is the approximation of the
c                               solution at the point (t,x).
c
        double precision        ux(npde)
c                               ux(1:npde) is the approximation of the
c                               spatial derivative of the solution at
c                               the point (t,x).
c
c output:
        double precision        dbdu(npde,npde)
c                               dbdu(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the unknown function u.
c
        double precision        dbdux(npde,npde)
c                               dbdux(i,j) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to the j-th component
c                               of the spatial derivative of the
c                               unknown function u.
c
        double precision        dbdt(npde)
c                               dbdt(i) is the partial derivative
c                               of the i-th component of the vector b
c                               with respect to time t.
c-----------------------------------------------------------------------
c loop indices:
        integer                 i, j
c-----------------------------------------------------------------------
c
c     assign dbdu(1:npde,1:npde), dbdu(1:npde,1:npde), and dbdt(1:npde)
c     according to the right boundary condition equation in terms of
c     u(1:npde), ux(1:npde), uxx(1:npde).
c
      do i = 1, npde
         do j = 1, npde
            dbdu(i,j) = 0.0d0
            dbdux(i,j) = 0.0d0
         end do
         dbdt(i) = 0.d0
		 dbdux(i,i) = 1.0d0
      end do
c
      return
      end
c-----------------------------------------------------------------------
      subroutine uinit(x, u, npde)
c-----------------------------------------------------------------------
c purpose:
c       this subroutine is used to return the npde-vector of initial
c       conditions of the unknown function at the initial time t = t0
c       at the spatial coordinate x.
c
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        double precision        x
c                               the spatial coordinate.
c
        integer                 npde
c                               the number of pdes in the system.
c
c output:
        double precision        u(npde)
c                               u(1:npde) is vector of initial values of
c                               the unknown function at t = t0 and the
c                               given value of x.
c-----------------------------------------------------------------------
      double precision pi
      parameter       (pi = 3.14159265358979323846264d0)
c-----------------------------------------------------------------------
      double precision a, b, w1, w2, gamma, xi, eta
	  common /tumor/ a, b, w1, w2, gamma, xi, eta
c-----------------------------------------------------------------------
c loop indices:
        integer                 i
c-----------------------------------------------------------------------
c
c     assign u(1:npde) the initial values of u(t0,x).
c
      do i = 1, npde
         u(i) = (exp(-((x - xi) * (x - xi))/(eta*eta)))/(eta*sqrt(pi))
      end do
c
      return
      end
c-----------------------------------------------------------------------
 
      subroutine header(nout)
c-----------------------------------------------------------------------
c purpose:
c       this subroutine writes a header describing the npde dimensional
c       parabolic partial differential equation
c                        ut = f(t, x, u, ux, uxx).
c-----------------------------------------------------------------------
      double precision a, b, w1, w2, gamma, xi, eta
	  common /tumor/ a, b, w1, w2, gamma, xi, eta
c-----------------------------------------------------------------------
c subroutine parameters:
c input:
        integer                 nout
c                               nout is the output unit number.
c-----------------------------------------------------------------------
c constants:
        double precision        t0
        parameter              (t0 = 0.0d0)
c
c-----------------------------------------------------------------------
c
      write(nout,95) 'Tumor growth equation:'
      write(nout,95) 'pde:'
      write(nout,95) '   u_t =  D*u_xx Dx*u_x,  '
      write(nout,95) 'domain:'
      write(nout,96) '   t0 =', t0, ' < t,'
      write(nout,96) '   xa =', a, ' <= x <= xb =', b, ','
c
      return
   95 format(a)
   96 format(a,e13.5,a,e13.5,a,e13.5,a,e13.5,a)
      end


