Attempting to perform the Runge-Kutta-Fehlberg Algorithm

I am attempting to code the Runge-Kutta-Fehlberg Algorithm for the solution of differential equations. I feel like my code is close, however on the example problem my code runs only 4 iterations and my updated h-step size exceeds the minimum value without having computed a low enough error estimation. Does anyone see the error in my code that would lead to perhaps R, h, or del being computed incorrectly? Perhaps my K's have something incorrect? I have stared at it for hours and can't locate it:

(I know my codes probably not following all the correct coding rules, but I'm really just looking to get the right solution and don't care about optimization apart from that)

syms t y(t)
eqn = diff(y,t) == (y^2+y)/t;
cond = y(1)==-2;
sol(t) = dsolve(eqn,cond)

sol(t) = 

%%%%%%%%%%%%%%%%%%%% 2c  %%%%%%%%%%%%%%%%%%%%%
a = 1;
b = 3;
j = 1;
TOL = 10.^(-4);
hmax = .5;
hmin = .02;
h = hmax;
FLAG = 1;
t = a:hmin:b;
y = zeros(size(t));
y_4 = zeros(size(t));
tempy_4 = zeros(size(t));
y_5 = zeros(size(t));
tempy_4(1) = -2;
y_4(1) = -2;
y_5(1) = -2;
n = numel(y);
i = 1;


a1 = 1 /4;
b1 = 3 /8;
b2 = 3 /32;
b3 = 9 /32;
c1 = 12 /13;
c2 = 1932 /2197;
c3 = 7200 /2197;
c4 = 7296 /2197;
d1 = 439 /216;
d2 = 8;
d3 = 3680 /513;
d4 = 845 /4104;
e1 = 1 /2;
e2 = 8 /27;
e3 = 2;
e4 = 3544 /2565;
e5 = 1859 /4104;
e6 = 11 /40;
while FLAG == 1
    
    fprintf('while statement begun \n');
    f = @(t,y) (y.^2 + y)./t;
    K1 = h *f(t(i),           tempy_4(i));
    K2 = h *f((t(i)+a1 *h),  (tempy_4(i) + a1 *K1));
    K3 = h *f((t(i)+b1 *h),  (tempy_4(i) + b2 *K1 + b3 *K2));
    K4 = h *f((t(i)+c1 *h),  (tempy_4(i) + c2 *K1 - c3 *K2 + c4 *K3));
    K5 = h *f((t(i)+h),      (tempy_4(i) + d1 *K1 - d2 *K2 + d3 *K3 - d4 *K4));
    K6 = h *f((t(i)+e1 *h),  (tempy_4(i) - e2 *K1 + e3 *K2 - e4 *K3 + e5 *K4 - e6 *K5));
    tempy_4(i+1) = y_4(i) + (25 /216) *K1 + (1408 /2565) *K3+(2197 /4104) *K4-(1 /5) *K5;
    y_5(i+1) = y_5(i)+(16 /135) *K1 + (6656 /12835) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;
    R = abs(y_5(i+1)-tempy_4(i+1)) /h;
    
    
    if R <= TOL
        fprintf('Step 2 \n');
        t(i+1) = t(i) + h;
        y_4(i+1) = y_4(i) + 25 /216 *K1 + 1408 /2565 *K3+2197 /4104 *K4-1 /5 *K5;
        i = i+1;
    end
    del = .84.*(TOL./R).^(.25);
    if del <= .1
        fprintf('Step 3.1 \n');
        h = .1*h;
    elseif del >= 4
        fprintf('Step 3.2 \n');
        h = 4.*h;
    else
        fprintf('Step 3.3 \n');
        h = del.*h;
    end
    if h > hmax
        fprintf('Step 4 \n');
        h = hmax;
    end
    if t(i) >= b
        fprintf('Step 5.1 \n');
        FLAG = 0;
    elseif t + h > b 
        fprintf('Step 5.2 \n');
        h = b - t;
    elseif h < hmin
        fprintf('minimum h exceeded \n');
        FLAG = 0;
    end
    j = j + 1;
    if j > 10
        FLAG = 0;
    end
end

while statement begun 

Step 3.3 

while statement begun 

Step 3.3 

while statement begun 

Step 3.3 

while statement begun 

Step 3.3 

minimum h exceeded 

plot(t,y_4)

Answer 1

[Torsten]

A partial answer is given by the code below - your test ODE is solved correctly with fixed time stepsize.

Errors in your original code that are corrected are:

h = hmin;

instead of

h = hmax;

y_5(i+1) = y_4(i)+(16 /135) *K1 + (6656 /12835) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;

instead of

y_5(i+1) = y_5(i)+(16 /135) *K1 + (6656 /12835) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;

12825 

instead of

12835

in the denominator (see answer by Jim Riggs)

The adaptive version gives wrong results. The main reason is that the time vector t is no longer valid because you change h from step to step. Thus evaluating your function with the t(i) arguments is false.

I suggest you start from the code below and try to incorporate the adaptive stepsize control step by step.

syms t y(t)
eqn = diff(y,t) == (y^2+y)/t;
cond = y(1)==-2;
sol(t) = dsolve(eqn,cond)

sol(t) = 

f_ana = matlabFunction(sol);
%%%%%%%%%%%%%%%%%%%% 2c  %%%%%%%%%%%%%%%%%%%%%
a = 1;
b = 3;
j = 1;
TOL = 10.^(-4);
hmax = .5;
hmin = .02;
h = hmin;
FLAG = 1;
t = a:hmin:b;
y = zeros(size(t));
y_4 = zeros(size(t));
tempy_4 = zeros(size(t));
y_5 = zeros(size(t));
tempy_4(1) = -2;
y_4(1) = -2;
y_5(1) = -2;
n = numel(y);
i = 1;

a1 = 1 /4;
b1 = 3 /8;
b2 = 3 /32;
b3 = 9 /32;
c1 = 12 /13;
c2 = 1932 /2197;
c3 = 7200 /2197;
c4 = 7296 /2197;
d1 = 439 /216;
d2 = 8;
d3 = 3680 /513;
d4 = 845 /4104;
e1 = 1 /2;
e2 = 8 /27;
e3 = 2;
e4 = 3544 /2565;
e5 = 1859 /4104;
e6 = 11 /40;
%while FLAG == 1

fprintf('while statement begun \n');

while statement begun 

f = @(t,y) (y.^2 + y)./t;
for i = 1:numel(t)-1
    K1 = h *f(t(i),           y_4(i));
    K2 = h *f((t(i)+a1 *h),  (y_4(i) + a1 *K1));
    K3 = h *f((t(i)+b1 *h),  (y_4(i) + b2 *K1 + b3 *K2));
    K4 = h *f((t(i)+c1 *h),  (y_4(i) + c2 *K1 - c3 *K2 + c4 *K3));
    K5 = h *f((t(i)+h),      (y_4(i) + d1 *K1 - d2 *K2 + d3 *K3 - d4 *K4));
    K6 = h *f((t(i)+e1 *h),  (y_4(i) - e2 *K1 + e3 *K2 - e4 *K3 + e5 *K4 - e6 *K5));
    y_4(i+1) = y_4(i) + (25 /216) *K1 + (1408 /2565) *K3+(2197 /4104) *K4-(1 /5) *K5;
    y_5(i+1) = y_4(i)+(16 /135) *K1 + (6656 /12825) *K3+(28561 /56430) *K4 - (9 /50) *K5+(2 /55) *K6;
    %R = abs(y_5(i+1)-tempy_4(i+1)) /h;
    
    
    %if R <= TOL
    %    fprintf('Step 2 \n');
    %    t(i+1) = t(i) + h;
    %    y_4(i+1) = y_4(i) + 25 /216 *K1 + 1408 /2565 *K3+2197 /4104 *K4-1 /5 *K5;
    %    i = i+1;
    %end
    %del = .84.*(TOL./R).^(.25);
    %if del <= .1
    %    fprintf('Step 3.1 \n');
    %    h = .1*h;
    %elseif del >= 4
    %    fprintf('Step 3.2 \n');
    %    h = 4.*h;
    %else
    %    fprintf('Step 3.3 \n');
    %    h = del.*h;
    %end
    %if h > hmax
    %    fprintf('Step 4 \n');
    %    h = hmax;
    %end
    %if t(i) >= b
    %    fprintf('Step 5.1 \n');
    %    FLAG = 0;
    %elseif t + h > b 
    %    fprintf('Step 5.2 \n');
    %    h = b - t;
    %elseif h < hmin
    %    fprintf('minimum h exceeded \n');
    %    FLAG = 0;
    %end
    %j = j + 1;
    %if j > 10
    %    FLAG = 0;
    %end
end
hold on
plot(t,y_4)
plot(t,y_5)
plot(t,f_ana(t))
hold off
grid on

Answer 2

[Jim Riggs]

I found the error indicated below.

Note that your notation makes your code much harder to verify because you have separated the coefficient magnitudes from their signs. For example, C3 has a magnitude of 7200/2197, and a negative sign. So I have to look at two different places to verify that C3 is implemented correctly. If all of the K equations were implemented with "+" for each term, then you put the sign on the coefficient, it is much easier to verify; C3 = -7200/2197.