% p25.m - stability regions for ODE formulas

% Runge-Kutta:
  clf
  z = exp(1i*pi*(0:200)/100); r = z-1;
  plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
  w = 0; W = w; for i = 2:length(z)                  % order 1
    w = w-(1+w-z(i)); W = [W; w]; end, plot(W)
  w = 0; W = w; for i = 2:length(z)                  % order 2
    w = w-(1+w+.5*w^2-z(i)^2)/(1+w); W = [W; w];
    end, plot(W)
  w = 0; W = w; for i = 2:length(z)                  % order 3
    w = w-(1+w+.5*w^2+w^3/6-z(i)^3)/(1+w+w^2/2); W = [W; w];
    end, plot(W)
  w = 0; W = w; for i = 2:length(z)                  % order 4
    w = w-(1+w+.5*w^2+w^3/6+w.^4/24-z(i)^4)/(1+w+w^2/2+w.^3/6);
    W = [W; w]; end, plot(W)
  axis([-5 2 -3.5 3.5]), axis square, grid on, title Runge-Kutta
