Hi,
I see that you have a 2-Step BDF and are looking to change it to a 4-Step BDF method. In order to do so, you need to adjust the coefficients for the BDF formula and provide three additional starting values (since the 4-step BDF method requires four previous values to compute the next value). These additional starting values can be obtained using a method with sufficient order, such as the 4th-order Runge-Kutta method.
This can be done by computing the first four values of 'y' using the 4th-order Runge-Kutta method followed by a 4-step BDF method to compute the subsequent values. The coefficients (48/25, -36/25, 16/25, -3/25) and (12/25) are to be used, specific to the 4-step BDF method.
Here is a generic code snippet along the lines of your existing code:
(Make sure to adjust the rest of your code to handle the output 'yend' appropriately)
%% 4-step BDF method
f = @(t, y) cos(2*t + y) + (3/2)*(t - y);
y0 = 1;
a = 0;
b = 1;
N0 = 10;
p = 4; % 4-step BDF
j = 1;
for k = 0:5
N = N0*2^k;
NN(j) = N;
h = (b-a)/N;
hh(j) = h;
tn = a:h:b;
yn = zeros(1,length(tn));
y(1) = y0;
%% Runge-Kutta for initial values
for i = 1:p-1
k1 = h*f(tn(i), y(i));
k2 = h*f(tn(i) + h/2, y(i) + k1/2);
k3 = h*f(tn(i) + h/2, y(i) + k2/2);
k4 = h*f(tn(i) + h, y(i) + k3);
y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
end
t = a + (p-1)*h;
for i = p:length(tn)-1
ynew = fzero(@(ynew) ynew - (...
(48/25)*y(i) - (36/25)*y(i-1) + (16/25)*y(i-2) - (3/25)*y(i-3) + ...
(12/25)*h*f(t+h, ynew)), y(i));
y(i+1) = ynew;
t = t + h;
end
plot(tn, y)
hold on;
yend(j) = y(end);
j = j+1;
end
figure(1);
xlabel('t');
ylabel('y(t)');
[mk, nk] = size(yend);
for k = 1:nk
% Add any code here if needed for post-processing
end
Hope this helps!
