% Usage: plotrk(r,a)
% Plot the boundary of the stability region of Runge-Kutta methods
%
% Input:
% r - Matlab inline function r(x,y)
% a - size of the region
%
% Output: Plot of the staility region in white
%
% Examples:
% plotrk(@myfunc,3);          here 'myfunc' is a user-defined function in M-file
% plotrk(inline('((1+x*0.7)^2+(y*0.7)^2)/((1-x*0.3)^2+(y*0.3)^2)','x','y'),
% 5);    theta-method with theta=0.7
% r=inline('((1+x*0.3)^2+(y*0.3)^2)/((1-x*0.7)^2+(y*0.7)^2)','x','y');
% theta-method with theta=0.7
% plotrk(r,3);
function []=plotrk(r,a)
[x,y] = meshgrid(-a:a/25:a,-a:a/25:a);
z=zeros(size(x));
for i=1:length(x)
    for j=1:length(x)
        z(i,j)=r(x(i,j),y(i,j));
    end
end
colormap copper;
contourf(x,y,z,[1,1]);