%Math 316 - Spring 2017
%Numerical ODEs
%MATLAB code for EXAMPLEs 5, 6, 7 
%Written by M. Motamed @ UNM 
%===============================================================
%Plotting stability regions of RK1, RK2, and RK4
%===============================================================
clear all

%RK1
theta=linspace(0,2*pi,1000);
figure
set(gca,'fontsize', 20); 
plot(1+cos(theta),sin(theta),'k-')
axis equal
hold on

%RK2
phi=linspace(-pi,pi,1000);
M=exp(1i*phi);
z1=1+sqrt(2*M-1);
z2=1-sqrt(2*M-1);
%figure
plot(z1,'k')
hold on
plot(z2,'k-')
axis([-1 3 -2 2])

%RK4
[x,y]=meshgrid(linspace(-4,4,1000),linspace(-4,4,1000));
z=x+1i*y;
%figure
contour(x,y,abs(1-z+(1/2)*z.^2-(1/6)*z.^3+(1/24)*z.^4),[1,1],'k-','linewidth',2)
axis([-2 4 -4 4])
grid
