%% SECTION 6.1 - EULER'S METHOD

%%
% Here I use the function *myeuler* (from pages 104-105 of _Differential 
% Equations with MATLAB_) implementing Euler's method to solve y' = 2y - 1. 
% It takes as input the function f, the initial time tinit, the initial 
% y-value yinit, the final time value b and the number of steps n. 

syms t y
f = @(t,y) 2*y - 1

%%
% First I use step size 0.1.

[t1,y1] = myeuler(f,0,1,0.4,4); [t1,y1]

%%
% Next I use step size 0.05.

[t2,y2] = myeuler(f,0,1,0.4,8); [t2,y2]

%%
% Finally I use step size 0.025

[t3,y3] = myeuler(f,0,1,0.4,16); [t3,y3]

%%
% The exact solution is 
y = dsolve('Dy = 2*y-1, y(0) = 1','t')

%%
% Here is a table of values.  The first entry is t, followed by the
% approximate solutions evaluated at that value of t, followed by the exact
% solution evaluated at that value of t.  I build the table T a row at a
% time using a *for* loop.  I initialize T with an empty row.

format short
T = [ ];
for k=0:4
    A = [t1(k+1),y1(k+1),y2(2*k+1),y3(4*k+1),subs(y,'t',t1(k+1))];
    T = [T;A];
end
T
%%
% I use *double* to convert the entries to decimals.
double(T)
%%
% I plot the three approximate solutions and the exact solution.  I do
% two types of plots, one where I show just the points computed by Euler's
% method, the other where the points are connected by line segments. 

plot(t1,y1,'.c','MarkerSize',15), hold on
plot(t2,y2,'.r','MarkerSize',15)
plot(t3,y3,'.g','MarkerSize',15)
ezplot(y,[0,0.4])
title 'Solutions using Eulers method and the exact solution'
legend1 = legend({'step size 0.1', 'step size 0.05', 'step size 0.025','exact'},...
    'Position',[0.1465 0.7205 0.2171 0.1869]);
hold off

%%

plot(t1,y1,'c'), hold on
plot(t2,y2,'r')
plot(t3,y3,'g')
ezplot(y,[0,0.4])
title 'Solutions using Eulers method and the exact solution'
legend1 = legend({'step size 0.1', 'step size 0.05', 'step size 0.025','exact'},...
    'Position',[0.1465 0.7205 0.2171 0.1869]);
hold off

%%
% I see from both the plots and the table that the error in the approximate 
% solutions decreases as the number of steps increases, but it is still 
% fairly large.

%%
% I can look at the errors. Here is a table with t in the first column,
% y1(t) in the second, y(t) in the third and e1(t)=y(t)-y1(t) in the
% fourth. Using *double* converts the entries in the table to decimals.

format short
T = [ ];
for k=0:4
    A = [t1(k+1),y1(k+1),subs(y,'t',t1(k+1)),subs(y,'t',t1(k+1))-y1(k+1)];
    T = [T;A];
end
double(T)
%%
% Here is a table with t in the first column,
% y3(t) in the second, y(t) in the third and e1(t)=y(t)-y3(t) in the
% fourth. I use the same t values as in the first table of errors.
format short
T = [ ];
for k=0:4
    A = [t3(4*k+1),y3(4*k+1),subs(y,'t',t3(4*k+1)),subs(y,'t',t3(4*k+1))-y3(4*k+1)];
    T = [T;A];
end
double(T)
%% 
% The errors are still increasing, but are much smaller than for y1.