

*********************************************************************** 
*     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 function which is being solved is y'(x) = f(x, y(x))
*     With f(x, y(x)) = const1*((x1)**(const2))*((x2)**(const3))
*     The user also inputs the initial 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.
*     2.0   Once i had a working system, i tightened up all the
*               loose ends, added in the ability for user controlled 
*               derivative function.  
*     3.0   Cleaned up the notation and variable names throughout.  
*               I also added in the beginnings of an adaptive step size
*
*     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.
***********************************************************************

***********************************************************************
*     MAIN PROGRAM
***********************************************************************

      program diffeq
      
      real*8 deltx, a,b,c, x, xtemp1, xtemp2, y, ytemp1, ytemp2, ytemp3
      real*8 xtemp3,deltx2,error
      integer j, nstep

      write(*,*) 'How many steps do you want to calculate? '
      read(*,*) nstep
      write(*,*) 'What increment do you want to use? '
      read(*,*) deltx
      write(*,*) ' What error limit do you want to use? '
      read(*,*) error
      write(*,*) 'What is the initial value of x, x0? '
      read(*,*) x0
      write(*,*) 'What is the initial value of y(x) at x = x0? '
      read(*,*) ytemp1
      write(*,*) 'What is the form of the function of x and y(x) '
      write(*,*) 'f(x,(y(x)) which equals the derivative of y(x)? '
      write(*,*) 'In the form dy(x)/dx = f(x,y(x)), with '
      write(*,*) ''
      write(*,*) 'f(x,y(x)) = ((a)*(x^b)*(y(x)^c)'
      read(*,*) a,b,c
      open(1, file='runge1.gnu')
      write(1,*) x0,ytemp1
      open(2, file='runge2.gnu')
      write(2,*) x0,ytemp1
***********************************************************************
*     Loop which creates file data points, calling on the sub runge
***********************************************************************

      ytemp2 = ytemp1

      do 60 j = 1, nstep
         xtemp1 = x0 + j*deltx

         call runge(xtemp1, ytemp1, deltx, a, b, c)

         write(1,*) xtemp1 , ytemp1           

 60   continue

      do 70 j = 1, nstep*2.0
         xtemp2 = x0 + j*deltx/2.0

         call runge(xtemp2, ytemp2, deltx/2.0, a, b, c)

         write(2,*) xtemp2 , ytemp2           

 70   continue

      close(1)
      stop
      end

***********************************************************************
*     Subroutine which calculates the value of y(x)
***********************************************************************

      subroutine runge(indvar, depvar, deltax, constant, xpower, ypower)
      
      real*8 deriv, h, indvar, deltax, depvar,k1, k2, k3, k4, constant
      real*8 xpower, ypower


      h = deltax / 2.0

      k1 = deltax * deriv(indvar,depvar,constant,xpower,ypower)
    
      k2 = deltax *
     *     deriv(indvar + h, depvar+k1 / 2.0, constant, xpower, ypower)
      
      k3 = deltax *
     *     deriv(indvar + h, depvar+k2 / 2.0, constant, xpower, ypower)
      
      k4 = deltax *
     *     deriv(indvar + deltax, depvar + k3, constant, xpower, ypower)
      
      depvar = depvar + (k1 + (2.0 * (k2 + k3)) + k4) / 6.0
      





      return
      end


***********************************************************************
*     Function which calculates the value of x2'(x1) = f(x1,x2(x1))
***********************************************************************

      function deriv(indvar,depvar,constant,xpower,ypower)
      
      real*8 deriv, indvar, depvar, constant, xpower, ypower


      deriv = constant*((indvar)**(xpower))*((depvar)**(ypower))

      return
      end
