echo on

% puzzler6, answered by Arend L. Schwab, TUDelft 2009
% Stability region for a Runge-Kutta 3 method
% There is a complete family of rk3 methods
% I used here the one used in ode23:
% k1 = f(tn,yn)
% k2 = f(tn+h/2,yn+h/2*k1)
% k3 = f(tn+3*h/4,yn+3*h/4*k2)
% yn+1 = yn + h/9*(2*k1+3*k2*4*k3)


x=-3:0.1:3;
y=-3:0.1:3;
[X,Y]=meshgrid(x,y);
L=X+Y.*i;
Z=abs(1+L+1/2.*L.^2+1/6.*L.^3);
contour(X,Y,Z,[1 1],'b')
xlabel('Re(h*\lambda)')
ylabel('Im(h*\lambda)')
title('Stability region for rk3')
axis equal
axis([-3 1 -3 3])
grid

% apparently you can numerically integrate a pure oscilatory system
% with this rk3 if h < 1.73/w0, where w0 is the highest eigenfrequency

echo off