rkck.f

Go to the documentation of this file.

       SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr)
       IMPLICIT NONE

       INTEGER n,NMAX
      double precision h,x,dydx(n),y(n),yerr(n),yout(n)
      parameter(nmax=50)
      INTEGER i
      double precision ak2(nmax),ak3(nmax),ak4(nmax),ak5(nmax),&
      ak6(nmax),ytemp(nmax),A2,A3,A4,A5,A6,B21,&
      B31,B32,B41,B42,B43,B51,B52,B53,&
      B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3,DC4,DC5,DC6
      parameter(a2=.2,a3=.3,a4=.6,a5=1.,a6=.875,b21=.2,b31=3./40.,&
      b32=9./40.,b41=.3,b42=-.9,b43=1.2,b51=-11./54.,b52=2.5,&
      b53=-70./27.,b54=35./27.,b61=1631./55296.,b62=175./512.,&
      b63=575./13824.,b64=44275./110592.,b65=253./4096.,c1=37./378.,&
      c3=250./621.,c4=125./594.,c6=512./1771.,dc1=c1-2825./27648.,&
      dc3=c3-18575./48384.,dc4=c4-13525./55296.,dc5=-277./14336.,&
      dc6=c6-.25)
      do i=1,n
        ytemp(i)=y(i)+b21*h*dydx(i)
      end do

      call source_population_eq(x+a2*h,ytemp,ak2)
      do i=1,n
        ytemp(i)=y(i)+h*(b31*dydx(i)+b32*ak2(i))
      end do
      call source_population_eq(x+a3*h,ytemp,ak3)
      do i=1,n
        ytemp(i)=y(i)+h*(b41*dydx(i)+b42*ak2(i)+b43*ak3(i))
      end do
      call source_population_eq(x+a4*h,ytemp,ak4)
      do i=1,n
        ytemp(i)=y(i)+h*(b51*dydx(i)+b52*ak2(i)+b53*ak3(i)+b54*ak4(i))
      end do
      call source_population_eq(x+a5*h,ytemp,ak5)
      do i=1,n
        ytemp(i)=y(i)+h*(b61*dydx(i)+b62*ak2(i)+b63*ak3(i)+b64*ak4(i)+&
      b65*ak5(i))
      end do
      call source_population_eq(x+a6*h,ytemp,ak6)
      do i=1,n
        yout(i)=y(i)+h*(c1*dydx(i)+c3*ak3(i)+c4*ak4(i)+c6*ak6(i))
      end do
      do i=1,n
        yerr(i)=h*(dc1*dydx(i)+dc3*ak3(i)+dc4*ak4(i)+dc5*ak5(i)+dc6*&
      ak6(i))
      end do
      return
      END