Here is a cleaned-up version of the Matlab script we developed in class on Monday implementing Euler’s method.

You should “step through” this code and make sure you understand what’s happening at each step (i.e., copy and paste the code line-by-line into the Matlab command window and examine what variables are created at each step).

As written below, the code does the computations for Example 3 in Section 2.7 of Boyce and DiPrima, i.e., for the initial value problem y’ = 4 – t – 2y, y(0) =1, with h = 0.1.  You change edit the code so it computes the Euler approximation for different h.  You should be able to recreate the results in Table 2.7.3. (For this, I would recommend creating saving the code below as an m-file, commenting out the line which set h, set h from the command line, and then run the script.)

Then you should try to change the code so it computes the Euler approximation for other initial value problems (such as Exercises 1, 3 or 11–again, try to recreate the numbers in the solutions).

A further exercise would be to plot the direction field for the differential equation on the same graph as the Euler approximation and exact solution.  Recall that Matlab code for producing direction fields can be found here.


%This script implements Euler's method
%for Example 2 in Sec 2.7 of Boyce & DiPrima


%For different differential equations y'=f(t,y), update in two places:
%(1) within for-loop for Euler approximations
%(2) the def'n of the function phi for exact solution (if you have it)


%also update step size h; initial conditions t0,y0; endpt t_end


%set parameters for Euler's method:
h = 0.1; %step size


%set initial conditions
t0 = 0;
y0 = 1;


%end point
t_end = 5;


%calculate number of steps
n = (t_end-t0)/h;


%t and y will be arrays containing Euler results
t(1) = t0;
y(1) = y0;


for k=1:n
yprime(k+1) = 4 - t(k) +2*y(k); %UPDATE: RHS is diff eqn for y'(t[k], y[k])
t(k+1) = t(k) + h;
y(k+1) = y(k) + yprime(k+1)*h;
end


%transpose t and y into column vectors
tvalues = t';
yvalues = y';


%UPDATE: create function phi for exact solution
phi = @(t) (-7/4) + 0.5*t + (11/4)*exp(2*t);


%plot Euler approximation and exact solution
plot(t,y);
hold on;
fplot(phi,[t0,t_end]);


%put values for t, Euler approximation, exact solution into single matrix
results = [tvalues yvalues phi(tvalues)]