*********************************************************************** 
*     Program diffeq.f
*     Author: Jesse Welch, for Jim Mahoney's Programming Workshop SP06
*     Date: 4/05/06
*     Version: 2.0
*
*     This program calculates the values of a function y(x) over the 
*     interval [x0,deltx*nstep], each of which the user is given full
*     choice.  The user also inputs the value of y(x) at x = x0.
*     In this way, the user can use this program to calculate any 
*     ordinary differential equation with initial values and then
*     using the outputted file, plot the values of y(x) over a range
*     supplied by the user.
*
*     Version History:
*     1.0 - Before the program was at all utilizable i went through
*               the stages of creating first the subroutine which would
*               carry out the calculation, followed by the deriv funct-
*               ion.  Once this was done, i followed up with a front
*               end so to speak with the user controls + file access.
*     
*     Bugs: None known of at this point, many testing cycles checked
*
*     Future: I am going to follow up with a batch file which runs and 
*               plots the outputted file. I also hope to extend my 
*               program to include an extension to systems of ordinary
*               differential equations linked to a partial differential
*               equation.  There is also a possibility to increase 
*               accuracy by including an adaptive step size sub.
***********************************************************************


      program diffeq
      
      real*8 x,y,deltx,const1,const2,const3
      integer i,j,nstep
      write(*,*) 'What do you want nstep to be? '
      read(*,*) nstep
      write(*,*) 'What do you want delta x to be? '
      read(*,*) deltx
      write(*,*) 'What is the value of y at t = 0? '
      read(*,*) y
      write(*,*) 'What is the value of the derivative of the '
      write(*,*) 'function in the derivative? '
      write(*,*) 'In the form dydx = f(x,y(x)), with '
      write(*,*) 'f(x,y(x)) = (a)*(x^b)*(y(x)^c)'
      read(*,*) const1, const2, const3
      open(1, file='runge.gnu')
      write(1,*) 0,y
      do 60 j = 1, nstep
         x = j*deltx
         call rk4(x,y,deltx,const1,const2,const3)
         write(1,*) x,y
 60   continue
      close(1)
      stop
      end

      subroutine rk4(x,y,deltx,const1,const2,const3)
      
      real*8 deriv, h, x, deltx, y
      real*8 k1, k2, k3, k4,const1,const2,const3
      integer i, N
      h= deltx/2.0
      k1 = deltx * deriv(x,y,const1,const2,const3)
      k2 = deltx * deriv(x+h,y+k1/2.0,const1,const2,const3)
      k3 = deltx * deriv(x+h,y+k2/2.0,const1,const2,const3)
      k4 = deltx * deriv(x+deltx, y+k3,const1,const2,const3)
      y = y + (k1+(2.*(k2 + k3))+ k4)/6.0
      return
      end

      function deriv(x1, x2,const1,const2,const3)
      
      real*8 deriv,x1,x2,const1,const2,const3
      integer i

      deriv=const1*((x1)**(const2))*((x2)**(const3))

      return
      end
