clear all
close all

echo on

% Problem Set D, Problem 1
% Part (a)
% Here is the correct way to enter Airy's equation.
airy = inline('[y(2); x*y(1)]', 'x', 'y');

[xfor, yfor] = ode45(airy, [0,2], [0,1]);
plot(xfor, yfor(:,1))
hold on
[xbak, ybak] = ode45(airy, [0,-2], [0,1]);
plot(xbak, ybak(:,1))

facA = dsolve('D2y = 0', 'y(0)=0', 'Dy(0)=1', 'x')
ezplot(facA, [-2, 2])
hold off
% The facsimile agrees well with the actual solution in a
% neighborhood of the origin.
pause
print -deps figd1-1

% Part (b)
% Now we can look at the solution for x near -16 = -4^2. 

[xbak, ybak] = ode45(airy, [0,-20], [0,1]);
plot(xbak, ybak(:,1), 'r')
hold on
axis([-18 -14 -1 1])

% After viewing the graph, we guess at initial conditions to 
% fit the curve more closely. We need to make this adjustment
% to match up with the apparent amplitude and frequency in this
% neighborhood.
pause

facB = dsolve('D2y = -16*y', 'y(-15.9)=0.6', 'Dy(-15.9)=0', 'x')
ezplot(facB, [-18 -14])
hold off
pause
print -deps figd1-2

% Part (c)
% Now we look at the solution for x near 16 = 4^2.

[xfor, yfor] = ode45(airy, [0,18], [0,1]);
plot(xfor, yfor(:,1), 'r')
hold on
axis([14 18 10^20 4*10^21])

% We know the solution should be close to a hyperbolic sine
% curve. We fit the parameters by estimation from the graph
% of the numerical solution.
pause

facC = '10^(-9)*sinh(4*x)'
ezplot(facC, [14 18])
hold off
pause 
print -deps figd1-3

% Part (d)

plot(xbak, ybak(:,1))
hold on
plot(xfor, yfor(:,1))
axis([-20 2 -2 3])
pause
print -deps figd1-4

% It appears that the frequency increases and the amplitude
% decreases as x goes to minus infinity. The increasing frequency
% is predicted by our analysis; the decreasing frequency is harder
% to explain.

echo off