subroutine ehpc(m,n,h,x,y,e,choice)
    implicit none

    ! Dummy variables
    integer :: m, n, choice
    double precision :: h, e
    double precision, dimension(1:n,0:m) :: y   ! Table of values
    double precision, dimension(0:m) :: x       ! Time t

    ! Other variables
    double precision, dimension(1:n) :: y0_i1   ! Initial estimate
    double precision, dimension(1:n) :: yc_i1   ! Modified initial estimate
    double precision, dimension(1:n) :: yip     ! Previous initial estimate
    double precision, dimension(1:n) :: yx_c    ! Current estimate (corrector)
    double precision, dimension(1:n) :: yx_n    ! Next estimate (corrector)
    double precision, dimension(1:n) :: et      ! Estimation of truncation error
    double precision, dimension(1:n) :: err     ! Current error
    double precision :: k1, k2, k3, k4, ek, s   ! Extra vars
    integer :: i, p

    yx_c = 0
    yx_n = 0
    yip  = 0
    do i = 4, m
        y0_i1 = 0
        ! Predictor equation

        ! Initial estimate
        do p=1, n
            call f(p, y(:, i-1), k1, choice)
            call f(p, y(:, i-2), k2, choice)
            call f(p, y(:, i-3), k3, choice)
            y0_i1(p) = y(p, i-4) + ((4/3)*h*(2*k1 - k2 + 2*k3))
        end do

        ! Improved initial estimate
        do p=1, n
            yc_i1(p) = y0_i1(p) + ((112/121)*(yip(p)))
        end do

        yx_c = yc_i1

        ! Corrector equation
        do
            ek=0
            s=0
            ! Compute current terms
            do p=1, n
                call f(p, yx_c,     k4, choice)
                call f(p, y(:,i-1), k1, choice)
                call f(p, y(:,i-2), k2, choice)
                yx_n(p) = ((9/8)*y(p,i-1)) - (y(p,i-3)/8) + ((3/8)*h*(k4+2*k1-k2))
            end do

            ! Compute error

            do p=1, n
                err(p) = abs((yx_n(p) - yx_c(p))/yx_n(p))
            end do

            ek = sum(err)/n
        !    print *, ek
            if (ek<e) then
                exit
            else
                yx_c = yx_n
                yx_n = 0
            end if
        end do

        ! Estimation of truncation error
        do p=1,n
            et(p) = (-9/121)*(yx_n(p) -y0_i1(p))
        end do

        ! Final estimate computation
        do p=1, n
            y(p, i) = yx_n(p) + et(p)
        end do
        x(i) = x(i-1) + h

        ! For next iteration
        do p=1, n
            yip(p) = yx_n(p) - y0_i1(p)
        end do
    end do
end subroutine ehpc