
%% Discussion of Stability 

%% default accuracy

dsolve('Dy = y + sin(2*x)', 'y(0)=-0.4', 'x')

f = inline('y+sin(2*x)','x','y');
[x,ya1]=ode45(f,[0,15],-0.4);
plot(x,ya1)
hold on
ezplot(' -0.2*sin(2*x)-0.4*cos(2*x) ',[0 15])
hold off
axis([0 15 -1 2])
xlabel 'x'
ylabel 'y'
title 'Fig. A. Exact and Approximate Solutions'

print -dpng figureA
%% increased accuracy

f = inline('y+sin(2*x)','x','y');
[x,ya1]=ode45(f,[0,20],-0.4);
plot(x,ya1)
hold on
options = odeset('AbsTol', 1e-8,'RelTol', 1e-5);
[x,ya2]=ode45(f,[0,20],-0.4, options);
plot(x,ya2)
ezplot(' -0.2*sin(2*x)-0.4*cos(2*x) ',[0 20])
hold off
axis([0 20 -1 2])
xlabel 'x'
ylabel 'y'
title 'Fig. B. Effect of Increased Accuracy'
print -dpng figureB

%% Example of a stable equation

g = inline('x.^2 - exp(y.^3)','x','y');
close all; axes; hold on
for y0= -3:0.5:1
   [x,y]=ode45(g,[0,3],y0);
   plot(x,y)
end
hold off
axis([0 3 -3 1.5])
xlabel 'x'
ylabel 'y'
title 'Fig. C. Numerical Solutions of a Stable Equation'

print -dpng figureC

