I need help to build Fortran 95 code that solves 3rd order ODE using Runge-Kutta-Fehlberg method. I am trying to solve y”’ = -2y”+y’+2y with y(0)=3, y'(0)=-2, and y”(0)=6. The exact solution should be

y=exp(-2t)+exp(-t)+exp(t)

The plot shows comparison of it with the computation.
My plot

And here is my code.

program ft

    !testing 45 rk method by examining y''' = -2y''+y'+2y, y''(0)=6, y'(0)=-2, and y(0)=3.
    !the exact solution should be y = exp(-2t)+exp(-t)+exp(t)
    implicit none
    integer, parameter :: i = SELECTED_REAL_KIND (10,200)
    REAL(i) :: t, val(3), dt, newval(3), newdt
    logical :: redo
    external runge

    t = 0.0

    !val = (y,y',y'')
    !initial condition
    val = (/3.0,-2.0,6.0/)
    !step size
    dt = 1.0E-4

    open (10, file = 'RK_test.dat', status = 'new')
    do while (t<4.0)
        redo = .true.
        write(10,*) t, val(1), val(2), val(3)

        do while (redo)
            call runge(t, dt, val, newdt, redo, newval)
            dt = newdt
        end do
        val = newval
        t = t + dt
    end do
    close(10)

end program ft

subroutine runge(t, dt, val, newdt, recalc, valp1)
    implicit none
    integer, parameter :: i = SELECTED_REAL_KIND (10,200)
    REAL(i), parameter :: deltacalc(2,6) = reshape([16.0/135, 0.0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55, &
                25.0/216, 0.0, 1408.0/2565,   2197.0/4104,  -1.0/5,  0.0], shape(deltacalc), order=[2,1])
    REAL(i), parameter :: coe(6,6) = reshape ([0.0,         0.0,          0.0,          0.0,         0.0,      0.0, &
                                           1.0/4,       0.0,          0.0,          0.0,         0.0,      0.0, &
                                           3.0/32,      9.0/32,       0.0,          0.0,         0.0,      0.0, &
                                           1932.0/2197, -7200.0/2197, 7296.0/2197,  0.0,         0.0,      0.0, &
                                           439.0/216,   -8.0,         3680.0/513,   -845.0/4104, 0.0,      0.0, &
                                           -8.0/27,     2.0,          -3544.0/2565, 1859.0/4104, -11.0/40, 0.0], &
                                           shape(coe), order=[2,1])
    REAL(i) :: valbar(3), k(6,3), delta(3), error(3), var(3), f(3)
    REAL(i), intent(in) :: t, dt, val(3)
    REAL(i), intent(out) :: newdt, valp1(3)
    logical, intent(out) :: recalc
    integer :: a, b, c, d
    External func

    k = 0
    valbar = val
    valp1 = val
    var = val

    do a = 1,6
        do b = 1,3
            do c = 1,6
                do d = 1,3
                    var(d) = var(d) + coe(a,c)*k(c,d)
                end do
            end do
            call func(var, f)
            k(a,b) = dt*f(b)
        end do
    end do
    recalc = .false.

    do b = 1, 3
        if (val(b)==0.0) then
            error(b) = 1.0E-3
        else
            error(b) = 10.0**(log10(abs(val(b)))-3.0)
        end if


        do a = 1, 6
            valbar(b) = valbar(b) + k(a,b)*deltacalc(1,a)
            valp1(b) = valp1(b) + k(a,b)*deltacalc(2,a)
        end do

        delta(b) = 0.84*(error(b)/(abs(valbar(b) - valp1(b))/dt))**(1.0/4) 
        
        if (abs(valbar(b) - valp1(b))/dt>error(b)) then
            recalc = .true.
            !print*, "oooops"
        end if
    end do


    newdt = minval(delta)*dt
    if (newdt > 0.01) then
        newdt = 0.01
    end if

end subroutine runge

subroutine func(var, out)
    implicit none
    integer, parameter :: i = SELECTED_REAL_KIND (10,200)
    REAL(i), intent(in) :: var(3)
    REAL(i), intent(out) :: out(3)

    out(1) = var(2) 
    out(2) = var(3)
    out(3) = -2*var(3) + var(2) + 2*var(1)
end subroutine func

This is linear ODE. Thus I expect that the computation match the solution much better. Please let me know if I am wrong about this.

Any suggestion would be appreciated.
Thank you in advance.