Matlab and the 3 Step Adams-Bashforth Method

Question

I am trying to run and plot the solutions to the 3-step Adams-Bashforth method and am unable to understand where my code is wrong. I am very new to Matlab and have been asked to code this without a good prior knowledge of Matlab. Matlab plots my exact solution fine on the interval but I am not having the same luck with my approximated solution. Below is my code and any help would be greatly appreciated.

% To solve y' = -3*y+6*t+5 s.t. y(0) = 2e^3-1 for -1 <= t <= 2 using 
% 3-step Adams-Bashforth method.

clear, clc, clf
f = @(t, y) -3*y+6*t+5; % RHS
h = 0.2; % time step size
T = 3; % This is length of the time interval for which you're solving for, i.e. 2-(-1) = 3.
N = T/h; % total number of times steps

F = @(t) 2*exp(-3*t)+2*t+1; % true solution

% Preallocations:
t = zeros(N+1, 1); 
y = zeros(N+1, 1);

% Initializations:
y(1) = 2*exp(3)-1; % Initial value of the ODE IVP

% Compute the solution at t = h by using the true solution:
t(2) = h;
y(2) = F(t(2));

t(3) = 2*h;
y(3) = F(t(3));

% Main loop for marching N steps:
for i = 2:N
t(i+2) = i*h; % time points
y(i+2) = y(i+1) + (h/12)*(23*f(t(i+1), y(i+1)) - 16*f(t(i), y(i)) +5*f(t(i-1),y(i-1))); % 2-step Adams-Bashforth method!!!
end

% Plotting:
plot(t, y), hold on % plot the numerical solution obtained by Adams-Bashforth

% Plot the true solution:
tt = linspace(-1, 2, 1000); % sampling points for true solution
ff = F(tt); % function values at sampling points
plot(tt, ff, 'r') % plot the true solution using the sampling values

legend('Adams-Bashforth 3-step', 'exact', 'location', 'nw') % adding legends

Below is the Graph I am currently outputting - which is very off for a 3rd order approximation.

Answer 1

Since that comment got a bit long here a "comment-answer"

• How does your solution look? Does it explode, does it converge to something different?
• Can you include a plot, if there is one where you can actually see something?
• As that is an explicit method: Does a reduction of the step-size help?
• y(0) = 2e^3-1 for -1 <= t <= 2 using What does this mean? Do you want to start at $0$ or at $-1$?

Edit:

If I understand correctly: You want to solve: \begin{align*} y'(t) &= -3y+6t+5 \qquad \text{in }[-1,2]\\ y_{-1} &= 2e^3-1 = y(-1) \end{align*}

But what your code actually solves is: \begin{align*} y'(t) &= -3y+6t+5 \qquad \text{in } [-1,2]\\ y_{0} &= 2e^3-1 \neq y(0) \end{align*}

Additionally:

Matlab will not let you have a negative value for y(t).

Yes, it is normal, that you can't enter negative numbers there. What you use in matlab is an array, and arrays start at index 1. Your solution is supposed to be on the intervall $[-1,2]$, so we have $$t_1 = -1,\qquad t_2 = -1+h \qquad t_3=-1+2h \qquad … \qquad t_i = -1+(i-1)h$$

In Matlab that would be:

 t[1] = -1 % Stores -1 at index 1
 t[2] = -1+h % Stores -1+h at index 2
 ...

So you should fix this line

t(i+2) = i*h; % time points

---

What you can try to test if the rest of your implementation is correct, is setting the initial value to: $$y_0=y(0) = 3$$ in line y(1) = 2*exp(3)-1; % Initial value of the ODE IVP

This would mean that you solve the correct ODE with the correct starting value at $t=0$.