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!!!

% 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.



1 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$?


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*}


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$.

  • $\begingroup$ Unsure as to how I can show the plot on mathstack, it follows the direction of the curve however is very volatile and doesn't follow the curve of the exact solution well. y(-1) = 2exp^3-1 is the solution given for the original value problem, however Matlab will not let you have a negative value for y(t). So I want to start at -1. I believe the error could lie where I get the 3rd value for y as I am unsure as to what t(3) and y(3) should be set equal to. $\endgroup$
    – user496388
    Oct 27, 2017 at 14:58
  • $\begingroup$ Simply take a screenshot from the plot-window. Then edit your question (there is an edit button below the question) and click on the image icon in the bar at the top of the edit window. There you can simply include your picture. $\endgroup$
    – P. Siehr
    Oct 27, 2017 at 15:01
  • $\begingroup$ @user496388 Have a look at the edit in my answer. $\endgroup$
    – P. Siehr
    Oct 27, 2017 at 15:20
  • $\begingroup$ I just changed my code and as I make the values of H smaller is works, thank you very much for your help! $\endgroup$
    – user496388
    Oct 27, 2017 at 15:33
  • $\begingroup$ Did you only change the size of h, or also what I mentioned in my answer? Since you are new: You can mark the answer that helped the most with a green arrow (on the left side of the answer). This will help other people with the same problem. Another indicator of good answers are up-votes. $\endgroup$
    – P. Siehr
    Oct 27, 2017 at 15:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.