%% Section 8.1 - Euler's Method

%% Why is it often necessary to use graphical and numerical techniques?

%%
% Here is an example of an ODE which MATLAB claims to solve explicitly, but the
% solution is very complicated. The initial value problem is
% 
% $$y'=y^2+t^2, \quad y(0)=1.$$
% 
% Here is MATLAB's "exact solution."
y = dsolve('Dy = y^2 + t^2, y(0)=1', 't'); pretty(y)
%%
% The solution is very complicated and involves $z$, which is the solution of
% a very complicated equation.   Also, it gives me an absurd plot:
ezplot(y)
%%
% There is a discussion of this example and how to solve it using MuPad on
% pp;. 72-74 of _Differential Equations with MATLAB_.
%%
% We can solve it exactly using Maple, with the command:
%
% dsolve({diff(y(t), t) = y(t)^2+t^2, y(0) = 1}, y(t))
%
% The solution is:
syms t
y = -t*(-(gamma(3/4)^2-pi)/gamma(3/4)^2*besselj(-3/4,1/2*t^2)+...
    bessely(-3/4,1/2*t^2))/(-(gamma(3/4)^2-pi)/gamma(3/4)^2*besselj(1/4,1/2*t^2)...
    +bessely(1/4,1/2*t^2));

%%
% To get a better idea of what the solution is like, we can view it
% graphically.  One method is to plot the solution.

ezplot(y,[0,0.8])
xlabel 't', ylabel 'y'
title 'Solution of dy/dt = y^2 + t^2'

%%
% Another method is to look at the direction field using the *quiver*
% command.  

[T,Y] = meshgrid(-3:0.2:3,-3:0.2:3);
S = T.^2 + Y.^2;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis tight
xlabel 't', ylabel 'y'
title 'Direction field for dy/dt = y^2 + t^2'

%%
% A third method is to solve numerically and plot the resulting numerical
% solution.  Here we do that using *ode45*.

f = @(t,y) y^2 + t^2;
ode45(f,[0,0.8],1)

%%
% Note that this plot looks the same as what we found when we plotted the
% explicit solution.

%%
% Here's another example where MATLAB finds a very complicated solution.
%%
% 
% $$y' = y^2 + t, \quad y(0)=0.$$
% 
y = dsolve('Dy = y^2 + t, y(0) = 0','t')

%%
% Here's the direction field.

[T,Y] = meshgrid(-3:0.2:3,-3:0.2:3);
S = T + Y.^2;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis tight
xlabel 't', ylabel 'y'
title 'Direction field for dy/dt = y^2 + t'

%%
% If we attempt to plot the exact solution, we run into a problem.  MATLAB
% doesn't seem to know what *airyAi* and *airyBi* are.  Here's a way around that.
% These functions come from MuPad and haven't been converted to MATLAB
% notation. (This example is discussed on pp. 61-63 of _Differential
% Equations with MATLAB_.)
tvals = 0:0.01:2;
plot(tvals,double(subs(y,'t',tvals)))
xlabel 't', ylabel 'y'
title 'Exact solution'

%%
% Can this be correct?

%%
% Alternatively we can plot the numerical solution obtained with *ode45*.

f = @(t,y) y^2 + t;
[t,ya] = ode45(f,[0,2],0); 
plot(t,ya)

%%
% Can this be correct?

%%
% Here is a problem MATLAB can't solve exactly.
%%
% 
% $$y' = \mbox{--}\,t^2 y^3 + \cos(t), \quad y(0) = 0.$$
% 
% There is no choice but to use graphical and numerical techniques to
% analyze this equation.  Here is the direction field.

[T,Y] = meshgrid(-3:0.2:3,-3:0.2:3);
S = -T.^2*Y.^3  + cos(T);
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis tight
xlabel 't', ylabel 'y'
title 'Direction field for dy/dt = -t^2 y^3 + cos(t)'

%%
% Here is a plot of the numerical solution obtained with *ode45*.

f = @(t,y) -t^2*y^3 + cos(t);
[t,ya] = ode45(f,[0:100],0);
plot(t,ya)

%%
% From this we see that the solution appears to tend slowly towards 0 in an
% oscillating manner as x gets very large.