% rkStability
% plots stability regions for order 1 through 4
%   explicit RK methods

%  solve F(w) = R(w) - z^p = 0
%  w = w - F(w)/F'(w)

     clear, home, close all

     z = exp(i*linspace(0,2*pi,400));

     hold on
   
     w = 0; W = w; 
     for i = 2:length(z)         % order 1
       w = w-(1+w-z(i)); 
       W = [W; w]; 
    end
    plot(W,'b')

    
     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,'g')
       
      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,'r')
      
     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,'k')
     axis equal


  xlabel('real(\lambda \Delta t)'); ylabel('imag(\lambda \Delta t)');