!Cash-Karp Runge Kutta
La teoria

module cashkarp
    implicit none
    integer,private,parameter::q=selected_real_kind(14,300)
    real(kind=q),private,parameter::n11776=11776,n16632=16632,n6237=6237,&
        n24948=24948,n1540=1540,n3262=3262,n37800=37800,n4600=4600, &
        n44275=44275,n6831=6831,n8855=8855,n9361=9361,n38500=38500,&
        n20125=20125,n27648=27648,n10240=10240,n39550=39550,&
        n148600=148600,n94675=94675,n7479=7479,n96768=96768,&
        n2530=2530,n10=10
    real(kind=q),private,parameter::a21=n11776*n16632
    real(kind=q),private,parameter::a31=n11776*n6237,a32=3*n11776*n6237
    real(kind=q),private,parameter::a41=n11776*n24948,a42=-3*n11776*n24948,&
        a43=4*n11776*n24948
    real(kind=q),private,parameter::a51=-11*n11776*n1540,a52=135*n11776*n1540,&
        a53=-140*n11776*n1540,a54=70*n11776*n1540
    real(kind=q),private,parameter::a61=n8855*n3262,a62=n8855*n37800,&
        a63=n8855*n4600,a64=n8855*n44275,a65=n8855*n6831
    real(kind=q),private,parameter::b51=n10240*n9361,b53=n10240*n38500,&
        b54=n10240*n20125,b56=n10240*n27648
    real(kind=q),private,parameter::b41=n2530*n39550,b43=n2530*n148600,&
        b44=n2530*n94675,b45=n2530*n7479,b46=n2530*n96768
    real(kind=q),private,parameter::b31=344565760,b33=-362700800,b34=997427200
    real(kind=q),private,parameter::b21=-1468938240,b22=n10*244823040
    real(kind=q),private,parameter::b11=979292160
    real(kind=q),private,parameter::c2=195858432,c3=293787648,c4=587575296
    real(kind=q),private,parameter::c5=979292160,c6=856880640
    !
    real(kind=q),dimension(:),private,allocatable::yp,ytn  !  corrisponde a k1.
    real(kind=q),dimension(:),private,allocatable::k2,k3,k4,k5,k6
    integer,private::nv=0
    private:: esempio_cashkarp
    public:: program_test,xcontrollo,cash_karp,hermite,xscashkarp
    !
contains
    subroutine xcontrollo(tot)
        !   Se il valore di tot e' zero non ci sono errori evidenti.
        real(kind=q),intent(out)::tot
        tot=abs(a21-c2)
        tot=tot+abs(a31+a32-c3)
        tot=tot+abs(a41+a42+a43-c4)
        tot=tot+abs(a51+a52+a53+a54-c5)
        tot=tot+abs(a61+a62+a63+a64+a65-c6)
        tot=tot+abs(b51+b53+b54+b56-b11)
        tot=tot+abs(b41+b43+b44+b45+b46-b11)
        tot=tot+abs(b31+b33+b34-b11)
        tot=tot+abs(b21+b22-b11)
        return
    end subroutine xcontrollo
    !
    subroutine xscashkarp(y,nm,not_ok)
        real(kind=q),intent(in),dimension(1:)::y
        integer,optional,intent(in)::nm
        logical,intent(out)::not_ok
        not_ok=.false.
        if(nv==0) then
            nv=size(y,1)
            if(present(nm)) then
                nv=max(nv,nm)
            end if
            allocate(ytn(nv),k2(nv),k3(nv),k4(nv))
            allocate(k5(nv),k6(nv),yp(nv))
            ytn(:)=0.0_q
            k2(:)=0.0_q
            k3(:)=0.0_q
            k4(:)=0.0_q
            k5(:)=0.0_q
            k6(:)=0.0_q
            yp(:)=0.0_q
        else if(size(y,1)>size(ytn,1)) then
            not_ok=.true.            
        end if
        return
    end subroutine xscashkarp
    !
    subroutine cash_karp(yn,nv,y,hv,f,yt)
        !
        !  nv   rappresenta il numero di variabili da integrare e deve
        !       valere almeno 2 perche' la prima variabile e' quella
        !       indipendente.
        !  hv   rappresenta l'ampiezza del passo di integrazione
        !  y(:) rappresenta il vettore delle variabili a inizio passo.
        !  yn(:) rappresenta il vettore delle variabili a fine passo.
        !  yt(:) rappresenta il vettore delle derivate prime a inizio passo
        !       ma e' opzionale perche', se non viene dato, viene 
        !       calcolato internamente. 
        !  ytn(:) rappresenta il vettore delle derivate prime a fine passo, che
        !       viene calcolato solo se yt(:) e' passato alla subroutine.
        !  f(z,w) e' la subroutine che calcola il vettore delle derivate prime
        !       quando  noto il vettore delle variabili.
        !  yp(:) rappresenta il vettore delle variabili a fine passo ma 
        !       calcolato con una precisione minore. Se yp(:) differisce
        !       eccessivamente da yn(:) bisogna rifare il passo riducendo
        !       il valore di hv. Si noti che yp non  in argomento poiche'
        !       viene data come memoria del modulo.
        !
        integer,intent(in)::nv
        real(kind=q),intent(in),dimension(1:)::y
        real(kind=q),intent(in),dimension(1:),optional::yt
        real(kind=q),intent(in)::hv
        real(kind=q),intent(out),dimension(1:)::yn
        interface
            subroutine f(z,w)
                integer,parameter::q=selected_real_kind(15,300)
                real(kind=q),intent(in),dimension(1:)::w
                real(kind=q),intent(out),dimension(1:)::z
            end subroutine f
        end interface
        real(kind=q)::h
        h=hv/b11
        if(present(yt)) then
            ytn(1:nv)=yt(1:nv)
        else
            call f(ytn,y)
        end if
        yp(1:nv)=y(1:nv)+(h*a21)*ytn(1:nv)
        call f(k2,yp)
        yp(1:nv)=y(1:nv)+(h*a31)*ytn(1:nv)+(h*a32)*k2(1:nv)
        call f(k3,yp)
        yp(1:nv)=y(1:nv)+(h*a41)*ytn(1:nv)+(h*a42)*k2(1:nv)+&
                 (h*a43)*k3(1:nv)
        call f(k4,yp)
        yp(1:nv)=y(1:nv)+(h*a51)*ytn(1:nv)+(h*a52)*k2(1:nv)+&
                 (h*a53)*k3(1:nv)+(h*a54)*k4(1:nv)
        call f(k5,yp)
        yp(1:nv)=y(1:nv)+(h*a61)*ytn(1:nv)+(h*a62)*k2(1:nv)+&
                 (h*a63)*k3(1:nv)+(h*a64)*k4(1:nv)+(h*a65)*k5(1:nv)
        call f(k6,yp)
        yn=y !   per trascrivere eventuali dati oltre i primi nv.
        yn(1:nv)=yn(1:nv)+(h*b51)*ytn(1:nv)+(h*b53)*k3(1:nv)+&
                 (h*b54)*k4(1:nv)+(h*b56)*k6(1:nv)
        !  Valore meno preciso usabile solo per controllo:
        yp(1:nv)=y(1:nv)+(h*b41)*ytn(1:nv)+(h*b43)*k3(1:nv)+&
                 (h*b44)*k4(1:nv)+(h*b45)*k5(1:nv)+(h*b46)*k6(1:nv)
        if(present(yt)) then
            call f(ytn,yn)
        end if       
        return
    end subroutine cash_karp
    !
    subroutine hermite(yora,y,yt,yn,ytn,theta,n)
        !
        !   interpolazione cubica (di Hermite) 
        !
        integer,intent(in)::n
        real(kind=q),intent(in)::theta
        real(kind=q),intent(in),dimension(1:)::y,yt,yn,ytn
        real(kind=q),intent(out),dimension(1:)::yora
        real(kind=q):: h,d1,d2,d3,d4,ds,dsq,tq
        !
        !   La variabile indipendente deve essere la prima.
        !
        h=yn(1)-y(1)
        tq=theta*theta
        ds=theta-1.0_q
        dsq=ds*ds
        d1=dsq*(2.0_q*theta+1.0_q)
        d2=dsq*theta*h
        d3=tq*(3.0_q-2.0_q*theta)
        d4=tq*ds*h
        yora(1:n)=d1*y(1:n)+d2*yt(1:n)+d3*yn(1:n)+d4*ytn(1:n)
        return
    end subroutine hermite
    !
    subroutine esempio_cashkarp(z,w)
        !
        !   Il tempo e' ottenuto risolvendo l'equazione differenziale:
        !   `t = 1 e il valore del tempo e' memorizzato nella prima
        !   componente di w. Il vettore z rappresenta il valore
        !   delle derivate prime corrispondenti ai valori delle
        !   variabili memorizzate nel vettore w.
        !
        real(kind=q),intent(in),dimension(1:)::w
        real(kind=q),intent(out),dimension(1:)::z
        real(kind=q)::r32
        !
        if(w(1)>=0.0_q) then
            z(1)=1                  ! `t 
            z(2)=w(4)               ! `x
            z(3)=w(5)               ! `y
            r32=w(2)*w(2)+w(3)*w(3)     
            r32=1.0_q/(r32*sqrt(r32))
            z(4)=-w(2)*r32          ! `vx
            z(5)=-w(3)*r32          ! `vy
        else
            !
            !  Per poter inizializzare il calcolo , w(1) deve essere 
            !  negativo e il vettore z(:) rappresenta allora il vettore 
            !  delle variabili e non il vettore delle loro derivate prime.
            !  Al tempo w(1).lt.0 il valore assoluto di w(1) indica
            !  l'eccentricita' che deve essere maggiore di 0 e inferiore
            !  ad 1.
            !  Non e' ammessa una eccentricita' esattamente nulla perche'
            !  allora w(1) non sarebbe negativo ma qualsiasi valore pur
            !  piccolissimo va bene...
            !
            r32=min(0.999999_q,abs(w(1)))
            z(1)=0.0_q                         !  t=0
            z(2)=1-r32                         !  1-e
            z(3)=0.0_q                         
            z(4)=0.0_q
            z(5)=sqrt((1.0_q+r32)/(1.0_q-r32)) !  sqrt((1+e)/(1-e))
            write(unit=*,fmt="(1a,4f18.10)") "     ",z(2:)
        end if
        return
        !
        !   Questo esempio e' il problema del moto di un pianeta attorno al sole
        !   che occupa l'origine delle coordinate. Il semiasse maggiore e' sempre
        !   uno, se il moto e' ellittico ossia se l'eccentricita', mai negativa,
        !   risulta essere inferiore ad 1. Se valesse esattamente 1 il moto
        !   sarebbe parabolico. Se e indica l'eccentricita', allora a t=w(1)=0 
        !   bisogna imporre queste condizioni iniziali: z(1)=0, z(2)=1-e, z(3)=0, 
        !   z(4)=0, z(5)=sqrt((1+e)/(1-e)) ossia il moto inizia al perielio. 
        !
    end subroutine esempio_cashkarp
    !
    subroutine program_test(ece,passie,girie,posizioni)
        real(kind=q),intent(in),optional::ece
        integer,intent(in),optional::passie,girie
        real(kind=q),dimension(0:,1:),intent(out),optional::posizioni
        integer::passi,giri,j,kp
        real(kind=q)::ec,hv,t0,s,x
        real(kind=q),dimension(5)::y,yn,yt,yora
        logical::boh
        !----
        write(unit=*,fmt=*)"Esegue alcune prove di funzionamento",&
           " di questo modulo (module cashkarp)"
        write(unit=*,fmt=*)" "
        call xcontrollo(x)
        write(unit=*,fmt=*)"1) Se vale 0.0 ",&
            "i coefficienti sono forse giusti : ",x
        if(present(ece)) then
            ec=ece
        else
            write(unit=*,fmt="(a)",advance="no") &
               " Assegna l'eccentricita' dell'orbita "
            read(unit=*,fmt=*) ec
            ec=max(1.0e-30_q,min(0.999999_q,ec))
        end if
        if(present(passie)) then
            passi=passie
        else
            write(unit=*,fmt="(a)",advance="no") &
               " Assegna il numero di passi"
            read(unit=*,fmt=*) passi
        end if
        if(present(girie)) then
            giri=girie
        else
            write(unit=*,fmt="(a)",advance="no") &
               " Assegna il numero di giri"
            read(unit=*,fmt=*) giri
        end if
        !
        if(present(posizioni)) then
            posizioni=0.0_q
        end if
        !
        call xscashkarp(y,5,boh)
        if(boh) then
            write(unit=*,fmt=*)"Fallisce l'inizializzazione di cashkarp: invia e stop!"
            read(unit=*,fmt=*)
            stop
        end if
        !   hv*passi deve corrispondere ad un dodicesimo di giro ossia π/6
        hv=2.0_q*atan(1.0_q)/(3.0_q*max(1,passi))
        yn(:)=0
        yn(1)=-ec
        call esempio_cashkarp(y,yn)
         !  Ogni j=passi dovrebbe aver fatto un dodicesimo di giro.
        do kp=1,giri*12-12
            do j=1,passi
                call cash_karp(yn,5,y,hv,esempio_cashkarp)
                y=yn
            end do
        end do
        do kp=1,12
            do j=1,passi
                call cash_karp(yn,5,y,hv,esempio_cashkarp)
                y=yn
            end do
            write(unit=*,fmt="(i5,4f18.10)") kp,y(2:5)
        end do
        write(unit=*,fmt=*)"Ora il tempo vale:      ",y(1)
        t0=giri*atan(1.0_q)*8.0_q
        write(unit=*,fmt=*)"Il tempo voluto sarebbe:",t0
        read(unit=*,fmt=*)
        !  
        if(.not.present(posizioni)) then
            return
        end if
        !
        !   Esegue questi calcoli solo se deve fare una tabella con 
        !   le posizioni distanziate un centesimo di tempo ovvero
        !   un centesimo di radiante, dato che fa un giro in 6.28318....
        !
        call esempio_cashkarp(yt,y)
        call cash_karp(yn,5,y,hv,esempio_cashkarp,yt) 
        do kp=0,630
            procedo:do
                if(yn(1)>t0) then
                    exit procedo
                end if 
                y=yn
                yt(1:5)=ytn(1:5)
                call cash_karp(yn,5,y,hv,esempio_cashkarp,yt) 
            end do procedo
            s=(t0-y(1))/(yn(1)-y(1))
            call hermite(yora,y,yt,yn,ytn,s,5)
            posizioni(kp,1)=yora(2)
            posizioni(kp,2)=yora(3)
            !  read(unit=*,fmt=*)
            t0=t0+0.01_q
        end do
        return
    end subroutine program_test
    !
end module cashkarp
!
!

Per verificare scommenta il program seguente:

!
!program xcashkarp
!    use cashkarp
!    implicit none
!    call program_test()
!    write(unit=*,fmt="(1x,1a)",advance="no") "Un invio per finire: "
!    read(unit=*,fmt=*)
!    stop
!end program xcashkarp
!