Adam-Moulton four-step derivation

clist = 
      Table[((-1)^k ) 
          Integrate[Binomial[-s, k], {s, -1, 0}], {k, 0, 
          4}]

clist[[1]] g4 + clist[[2]] (g4 - g3) + 
      clist[[3]] (g4 - 2  g3 + g2) + 
      clist[[4]] (g4 - 3  g3 + 3  g2 - g1) + 
      clist[[5]] (g4 - 4  g3 + 6  g2 - 4  g1 + g0)

Expand[%]

Together[%]

Predictor-Corrector

First, rewrite the initialization code in a style that is easier to read. Here's what we had before.

f[t_, y_] :=  -y;
t0 = 0. ;
y0 = 1. ;
tf = 1;
n = 10;
h = (tf - t0)/n;
tylist = {{t0, y0}};

      
Do[
tylist = 
          Append[tylist, 
{t + 
                h, 
tylist[[-1, 2]] + 
                h*f[t + 0.5*h, 
                    tylist[[-1, 2]] + 
                      0.5*h*f[t, 
                          tylist[[-1, 
                              2]]]]
}], 
{t, t0, t0 + 2  h, h}];
tylist

Here's a cleaner rewrite.

f[t_, y_] :=  -y;
t0 = 0. ;
y0 = 1;
tf = 1;
n = 10;
h = (tf - t0)/n;
tylist = {{t0, y0}};
t = t0;
y = y0;

      
Do[
 
        y = y + h*f[t + 0.5*h, y + 0.5*h*f[t, y]]; 
t = 
          t + h; 
tylist = 
          Append[tylist, {t, y}], 
{i, 0, 
          2}];
tylist

Now expand this to a predictor-corrector method using Adam-Bashforth four-step for prediction and Adam-Moulton four-step for correction.

Clear[f, t0, y0, tf, n, h, tylist, t, t1, t2, t3, t4, y1, y2, y3, y4, 
        y4p];
f[t_, y_] :=  -y;
t0 = 0. ;
y0 = 1;
tf = 1;
n = 40;
h = (tf - t0)/n;
tylist = {{t0, y0}};
t = t0;
y = y0;
  
Do[
 
        y = y + h*f[t + 0.5*h, y + 0.5*h*f[t, y]]; 
t = 
          t + h; 
tylist = 
          Append[tylist, {t, y}], 
{i, 0, 
          2}];
  
t1 = t0 + h;
t2 = t1 + h;
t3 = t2 + h;
y1 = tylist[[-3, 2]];
y2 = tylist[[-2, 2]];
y3 = tylist[[-1, 2]];
  
Do[
  
y4p = 
          y3 + h/24*(55  f[t3, y3] - 59  f[t2, y2] + 37  f[t1, y1] - 
                  9  f[t0, y0]); 
  
t4 = t3 + h; 
y4 = 
          y3 + h/720*(251  f[t4, y4p] + 646  f[t3, y3] - 264  f[t2, y2] + 
                  106  f[t1, y1] - 
                  19  f[t0, y0]); 
tylist = 
          Append[tylist, {t4, y4}]; 
t0 = 
          t1; 
t1 = t2; 
t2 = 
          t3; 
t3 = t4; 
y0 = 
          y1; 
y1 = y2; 
y2 = 
          y3; 
y3 = y4, 
{i, 3, 
          n - 1}];
tylist // TableForm
ListPlot[tylist]

A second example.  For comparison, we first plot the exact solution.  We use Mathematica's DSolve to get an exact solution.  (Note this is working with the initial condition.  I'm not sure what went wrong when we did this in class.)

Clear[y, t]

exactsoln = 
      DSolve[{y'[t] == 10  y[t] (1 - y[t]), 
          y[0] == 1/10}, y[t], t]

Note that the result is a replacement rule so we use the following syntaxt to make a plot.

exactsolnplot = 
      Plot[y[t] /. exactsoln, {t, 0, 2}, PlotRange -> All]

Clear[f, t0, y0, tf, n, h, tylist, t, t1, t2, t3, t4, y1, y2, y3, y4, 
        y4p];
f[t_, y_] :=  10  y (1 - y);
t0 = 0. ;
y0 = 0.1;
tf = 2;
n = 20;
h = (tf - t0)/n;
tylist = {{t0, y0}};
t = t0;
y = y0;

      
Do[
 
        y = y + h*f[t + 0.5*h, y + 0.5*h*f[t, y]]; 
t = 
          t + h; 
tylist = 
          Append[tylist, {t, y}], 
{i, 0, 
          2}];

      
t1 = t0 + h;
t2 = t1 + h;
t3 = t2 + h;
y1 = tylist[[-3, 2]];
y2 = tylist[[-2, 2]];
y3 = tylist[[-1, 2]];

      
Do[
  
y4p = 
          y3 + h/24*(55  f[t3, y3] - 59  f[t2, y2] + 37  f[t1, y1] - 
                  9  f[t0, y0]); 
  
t4 = t3 + h; 
y4 = 
          y3 + h/720*(251  f[t4, y4p] + 646  f[t3, y3] - 264  f[t2, y2] + 
                  106  f[t1, y1] - 
                  19  f[t0, y0]); 
tylist = 
          Append[tylist, {t4, y4}]; 
t0 = 
          t1; 
t1 = t2; 
t2 = 
          t3; 
t3 = t4; 
y0 = 
          y1; 
y1 = y2; 
y2 = 
          y3; 
y3 = y4, 
{i, 3, 
          n - 1}];
tylist // TableForm
approxsolnplot = Show[exactsolnplot, ListPlot[tylist]]

For comparison, let's look at the Midpoint approximation.

Clear[f, t0, y0, tf, n, h, tylist, t, t1, t2, t3, t4, y1, y2, y3, y4, 
        y4p];
f[t_, y_] :=  10  y (1 - y);
t0 = 0. ;
y0 = 0.1;
tf = 2;
n = 40;
h = (tf - t0)/n;
tylist = {{t0, y0}};
t = t0;
y = y0;
Do[
 
        y = y + h*f[t + 0.5*h, y + 0.5*h*f[t, y]]; 
t = 
          t + h; 
tylist = 
          Append[tylist, {t, y}], 
{i, 0, 
          n - 1}];
tylist // TableForm
approxsolnplot = Show[exactsolnplot, ListPlot[tylist]]