Index Exceeds Array Bounds

Hello! I am writing a code for a question using adaptive RK method. I encountered a problem in the while loop, where the system shows "Index Exceeds Array Bounds". I do not understand where have I gone wrong. Appreciate your help!

Note: dy_dx is a function file I wrote separately

syms dydx y x
dydx = -0.6*y + 10*(exp((-(x-2)^2)/(2*(0.075^2)))) ;
disp('The ODE for Q3 is:')
disp(dydx)
h = 0.5;
e = 0.00005;
% y(0) = 0.5

n = 0; i = 1;
x = zeros(i+1);
y = zeros(i+1);

while x(i) < 4
    % y1 is the est y for full step size
    k1 = zeros(i);
    k1(i) = dy_dx(x(i),y(i));
    k2 = zeros(i);
    k2(i) = dy_dx((x(i) + (h/2)),(y(i) + (k1(i)*h/2)));
    k3 = zeros(i);
    k3(i) = dy_dx((x(i) + (h/2)),(y(i) + (k2(i)*h/2)));
    k4 = zeros(i);
    k4(i) = dy_dx((x(i) + h),(y(i) + (k3(i)*h)));
    
    y1 = zeros(i);
    y1(i) = y(i) + ((h/6)*(k1(i) + 2*k2(i) + 2*k3(i) + k4(i)));
    
    % y2 is the est y for half step size
    h_half = h/2;
    k_1 = zeros(i);
    k_1(i) = dy_dx(x(i),y(i));
    k_2 = zeros(i);
    k_2(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_1(i)*h_half/2)));
    k_3 = zeros(i);
    k_3(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_2(i)*h_half/2)));
    k_4 = zeros(i);
    k_4(i) = dy_dx((x(i) + h_half),(y(i) + (k_3(i)*h_half)));
    
    y2_1 = zeros(i);
    y2_1(i) = y(i) + ((h_half/6)*(k_1(i) + 2*k_2(i) + 2*k_3(i) + k_4(i)));
    y(i) = y2_1(i);
    x(i) = x(i) + h_half;
    
    k_a = zeros(i);
    k_a(i) = dy_dx(x(i),y(i));
    k_b = zeros(i);
    k_b(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_a(i)*h_half/2)));
    k_c = zeros(i);
    k_c(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_b(i)*h_half/2)));
    k_d = zeros(i);
    k_d(i) = dy_dx((x(i) + h_half),(y(i) + (k_c(i)*h_half)));
    
    y2 = zeros(i);
    y2(i) = y(i) + ((h_half/6)*(k_a(i) + 2*k_b(i) + 2*k_c(i) + k_d(i)));
    
    del = y2(i) - y1(i);
    E_a = del/15;
    y2_new = zeros(i);
    y2_new(i) = y2(i) + E_a;
    err = abs(((y2_new(i) - y1(i))/y2_new(i))*100);
    if err < e
        break
    end
    y = y2_new(i);
    x = x + h_half;
    n = n + 1;
    i = i + 1;    
end

disp('The estimated y value is:')
disp(y2_new(i-1))
disp('The estimated x value is:')
disp(x(i-1))
disp('The number of iteration is:')
disp(n)

plot(x,y2_new)
hold on
plot(x,y);
legend('Half Step Size', 'Full Step Size')
xlabel('x'); ylabel('y');
title('Visual Representation of the Iteration')

Answer 1

[James Tursa]

Normally I would have expected to see code that sets the next values of x and y, but I don't see it. I.e., I am looking for lines like the following

x(i+1) = x(i) + stuff
y(i+1) = y(i) + stuff

So that when you increment i

i = i + 1;

the next values of x(i) and y(i) are available.

The way it looks to me is you increment i at the bottom of the while loop and then immediately test x(i) < 4. But the x(i) element will at some point not exist because it didn't get created in the loop.

The fix: Pre-allocate x and y properly, and set x(i+1) and y(i+1) inside your loop before incrementing i.

[LEW JUN KIT]

Thank you very much, Mr. James! It's all ok now.