solve a second-order differential equation using finite differences, I keep encountering the 'Index exceeds the number of array

"When attempting to run my MATLAB code to solve a second-order differential equation using finite differences, I keep encountering the 'Index exceeds the number of array elements' error. I've tried adjusting the length of the Vb vector in the findiff.m file as suggested earlier, but the issue still persists. Are there any other suggestions I could try to resolve this problem?"

the error is

result =

1

Index exceeds the number of array elements. Index must

not exceed 1.

Error in findiff (line 19)

Vb(N - 1) = Vb(N - 1) + (1 - h / 2 * feval(p, Vt(N - 1))) * beta;

Error in main (line 9)

result = findiff(@p, @q, @r, a, b, alpha, beta, N);

main.m

% main.m

a = 0; % Batas bawah interval

b = 4; % Batas atas interval

alpha = 1.25; % Kondisi awal x(0)

beta = -0.95; % Kondisi awal x(4)

N = 20; % Resolusi (banyaknya langkah diskretisasi)

% Panggil fungsi findiff dengan menggunakan fungsi p, q, dan r dari file terpisah

result = findiff(@p, @q, @r, a, b, alpha, beta, N);

% Tampilkan hasil atau lakukan apa yang Anda inginkan dengan solusi

disp(result);

trisys.m

function X = trisys(A, D, C, B)

N = length(B);

% Eliminasi maju (forward elimination)

for k = 2:N

mult = A(k-1) / D(k-1);

D(k) = D(k) - mult * C(k-1);

B(k) = B(k) - mult * B(k-1);

end

% Substitusi mundur (back substitution)

X(N) = B(N) / D(N);

for k = N-1:-1:1

X(k) = (B(k) - C(k) * X(k+1)) / D(k);

end

end

findiff.m

function F = findiff(p, q, r, a, b, alpha, beta, N)

% Inisialisasi vektor T dan X

T = zeros(1, N + 1);

X = zeros(1, N - 1);

% Inisialisasi vektor Va, Vb, Vc, Vd, dan h

Va = zeros(1, N - 2);

Vb = zeros(1, N - 1);

Vc = zeros(1, N - 2);

Vd = zeros(1, N - 1);

h = (b - a) / N;

% Menghitung vektor titik-titik grid Vt

Vt = a+h:h:a+h*(N-1);

% Menghitung vektor Vb sesuai dengan persamaan

Vb = -h^2 * feval(r, Vt);

Vb(1) = Vb(1) + (1 + h / 2 * feval(p, Vt(1))) * alpha;

Vb(N - 1) = Vb(N - 1) + (1 - h / 2 * feval(p, Vt(N - 1))) * beta;

% Menghitung vektor Vd sesuai dengan persamaan

Vd = 2 + h^2 * feval(q, Vt);

% Menghitung vektor Va sesuai dengan persamaan

Vta = Vt(2:N - 1);

Va = -1 - h / 2 * feval(p, Vta);

% Menghitung vektor titik-titik grid Vtc

Vtc = Vt(1:N - 2);

% Menyelesaikan sistem persamaan tridiagonal dengan trisys

X = trisys(Va, Vd, Vc, Vb);

% Menggabungkan vektor T dan X untuk membentuk matriks F

T = [a, Vt, b];

X = [alpha, X, beta];

F = [T' X'];

end

p.m

% p.m

function result = p(t)

result = 2 * t / (1 + t^2);

end

q.m

% q.m

function result = q(t)

result = -2 / (1 + t^2);

end

r.m

% r.m

function result = r(t)

result = 1

end

Answer 1

[Torsten]

Replace your functions by

p.m

% p.m
function result = p(t)
    result = 2 * t ./ (1 + t.^2);
end

q.m

% q.m
function result = q(t)
    result = -2 ./ (1 + t.^2);
end

r.m

% r.m
function result = r(t)
    result = ones(size(t));
end