implicit none
      integer i,n
      real a,b,h,x,y,z,xl,yl,zl,f,g
      real yp,y0,zp,z0
      real k1,k2,k3,k4
      real l1,l2,l3,l4
      real exacty,exactz


C     Input parameters
C     n=number of steps
C     x(0)=a=starting point
C     x(n)=b=end point
C     y(0)=y0=initial value of y
C     z(0)=z0=initial value of z

      print*,' Enter n'
      read(5,*) n
      print*,' Enter x(0),x(n),y(0),z(0)'
      read(5,*) a,b,y0,z0

C     Calculate h = step-length
      h=(b-a)/float(n)


 100  continue
      y=y0
      z=z0
      x=a

      write(6,30) x,y,z
 30   format(/,' Solving with 4th Runge-Kutta gives:',/,
     &         '   X          Y         Z',/,
     &         1x,f5.3,1x,f10.6,1x,f10.6)


      do i=1,n
         xl=x
         yl=y
	 zl=z
         k1=h*f(xl,yl,zl)
         l1=h*g(xl,yl,zl)
         k2=h*f(xl+0.5*h,yl+0.5*k1,zl+0.5*l1)
         l2=h*g(xl+0.5*h,yl+0.5*k1,zl+0.5*l1)
         k3=h*f(xl+0.5*h,yl+0.5*k2,zl+0.5*l2)
         l3=h*g(xl+0.5*h,yl+0.5*k2,zl+0.5*l2)
         k4=h*f(xl+h,yl+k3,zl+l3)
         l4=h*g(xl+h,yl+k3,zl+l3)
         x=xl+h
         y=yl+(k1+2*k2+2*k3+k4)/6.0
         z=zl+(l1+2*l2+2*l3+l4)/6.0
         exacty=sin(x)
	 exactz=cos(x)
         write(6,2) x,y,z,exacty,exactz
      enddo

 2    format(1x,f5.3,1x,f10.6,1x,f10.6,1x,f10.6,1x,f10.6)

      stop
      end


C     Input the function f(x,y,z) which is dy/dx

      real function f(x,y,z)
      implicit none
      real x,y,z
      f=z
      return
      end


C     Input the function g(x,y,z) which is dz/dx

      real function g(x,y,z)
      implicit none
      real x,y,z
      g=-sin(x)
      return
      end

C     Here we are solving y''=-sin x
C     We have set z=y'   and z'=-sin x
C
C     The exact solution is y=sin x
C                           z=cos x
C     Be sure you start with x=a
C                            y=sin(a)
C                            z=cos(a)
C     e.g start with x=0, y=0, z=1. Perhaps integrate over pi with 20 steps.
C
C     A sample run is:
C  Enter n
C  20
C  Enter x(0),x(n),y(0),z(0)
C  0 3.14159 0 1
C
C Solving with 4th Runge-Kutta gives:
C   X          Y         Z
C 0.000   0.000000   1.000000
C 0.157   0.156434   0.987688   0.156434   0.987688
C 0.314   0.309016   0.951057   0.309017   0.951057
C 0.471   0.453990   0.891007   0.453990   0.891007
C 0.628   0.587784   0.809017   0.587785   0.809017
C 0.785   0.707106   0.707107   0.707106   0.707107
C 0.942   0.809016   0.587786   0.809017   0.587786
C 1.100   0.891005   0.453991   0.891006   0.453991
C 1.257   0.951055   0.309018   0.951056   0.309018
C 1.414   0.987687   0.156436   0.987688   0.156436
C 1.571   0.999999   0.000001   1.000000   0.000001
C 1.728   0.987688  -0.156433   0.987689  -0.156433
C 1.885   0.951056  -0.309016   0.951057  -0.309015
C 2.042   0.891006  -0.453989   0.891007  -0.453989
C 2.199   0.809017  -0.587784   0.809018  -0.587784
C 2.356   0.707107  -0.707106   0.707108  -0.707105
C 2.513   0.587786  -0.809016   0.587787  -0.809016
C 2.670   0.453992  -0.891006   0.453993  -0.891005
C 2.827   0.309018  -0.951056   0.309020  -0.951056
C 2.985   0.156436  -0.987689   0.156437  -0.987688
C 3.142   0.000002  -1.000001   0.000003  -1.000000