First download the codes and save them in files 
(with the same names) in the same directory, 
say in a directory called 
      NUMERICS

Then open MATLAB in the directory NUMERICS 
and type the commands (where >> is the MATLAB prompt) 
 
>> diary
>> format long

Recall that, after you type a command in MATLAB, you press "return" 
to execute it.  If you type a semicolon (;) after a command, 
the computer will execute it, but will not display the output.  


    EULER'S METHOD

Suppose that you want to solve the initial value problem (IVP) 

	dy
	-- = y - t^2 + 1 ,        0 <= t <= 2
	dt 

	y(0) = 0.5 ,

which is solved in Example 1 of Section 5.2 (page 259) 
and in Example 1 of Section 5.3 (page 268).  

To solve this IVP using Euler's method, 
you will use the function 
    euler(f,a,b,ya,N)
which is saved in the file  euler.m ; 
note that in MATLAB the name of the function 
and the name of the file must be the same 
(except, of course the standard MATLAB extension .m ).  
Here  f  is the name of the function  f(t,y) 
giving the right-hand side of the differential equation 
(in this case,  f(t,y) = y - t^2 + 1 ), 
a  and  b  are the initial and the final "times", 
ya  is the initial value  y(a)  of the function y, 
and  N  is the number of intervals in which 
you divide the interval  [a,b].  

The program  euler.m  does the following (look at the file).  
First it computes the stepsize  h 
- note that at the end of the line there is a semicolon (;) 
which tells MATLAB not to print on the screen 
the result of executing this line.  
Then the code defines the vector-rows  T  and  Y , 
each of length  N+1 , and initializes them with zeros.  
If you type 

>> x = zeros(1,3)

and press "return", MATLAB will type 

x =

     0     0     0


The command that follows in the code  euler.m , 

T=a:h:b;

initializes the mesh points as follows:  
    T(1)   = a
    T(2)   = a + h
    T(3)   = a + 2h
    ...
    T(N)   = a + (N-1)h
    T(N+1) = a + Nh = b
(note that in MATLAB the counter starts from  1 , 
so that the indices are shifted by one in comparison 
with the textbook).  

The next line in  euler.m  defines the initial value 
of the unknown function:

Y(1)=ya;

The loop below computes the values of the approximations 
to the solution at the mesh points:  

for j=1:N
   Y(j+1)=Y(j)+h*feval(f,T(j),Y(j));
end

The command 
    feval(f,T(j),Y(j));
returns the value of the function  f(t,y)  
for  t = T(j)  and  y = Y(j) .

Finally, the line

E=[T' Y'];

defines the value E that the program returns.  

To run  euler.m , you have to save the right-hand side 
of the differential equation, namely, 
     y - t^2 + 1 
in a file.  This function is saved in the file  rhs.m  
- again, the name of the file is the same as the value 
of the function it defines.  

To solve the IVP above by using Euler's method 
with  N = 10 , you have to type

>> euler(@rhs,0,2,0.5,10)

and press "return" - you will see the values 
of the approximations to the true solutions 
- compare these values with the ones 
from the table on page 259.  

The symbol  @  in front of the function's name 
is called the "handle" to the function rhs 
- this way you make MATLAB call the function 
rhs  within the routine  euler.m  every time 
it sees  f  in this code.  

Alternatively, you can type 

>> A = euler(@rhs,0,2,0.5,10);

- this will assign the new variable A 
the output of the program (this is what the 
variable  E  on the last line of the program is for).  
Then type 

>> A

and see what happens.  


    TAYLOR'S METHOD OF ORDER 2

To solve the IVP above by using Taylor's method of order 2, 
you have to run the MATLAB routine  taylor2  saved in the 
file  taylor2.m .

The function  taylor2  takes the following arguments:
     taylor2(df,a,b,ya,N)
- here df is a vector-row consisting of two elements, 
     df=[y' y'']
where  y'  is the derivative of the unknown function  y(t), 
which, of course, is equal to the right-hand side  f(t,y) 
of the differential equation, 
and  y''  is the second derivative of  y(t) , 
which is equal to total derivative with respect to  t 
of  f(t,y)  - please do not forget that y depends on  t !  

From Example 1 in Section 5.3 (page 268), 
you know that 
    y'  = f(t,y) = y-t^2+1 , 
and 
    y'' = (d/dt) f(t,y(t)) = y-t^2+1-2t .

Since the code  taylor2.m  needs both  y'  and  y'' , 
we will save them both in the file  twoders.m 
which, as you can guess, contains the function  twoders(t,y) .  
Every time this program is called, it returns two numbers 
- the first one equal to  y' , i.e., to  y-t^2+1 , 
and the second one equal to  y'' , i.e., to  y-t^2+1-2t .
In the MATLAB file  twoders.m  this is done as follows:

twoD = [y-t^2+1 y-t^2+1-2*t];

- note that there is a space between  y'  and  y'' , 
but there are no spaces in  y-t^2+1  and in  y-t^2+1-2*t 
- PLEASE don't forget this when you create the file 
threeders.m  which will return the derivatives  
y' ,  y'' , and  y'''  of the appropriate function.  

When you run  taylor2.m  as follows:

>> A = taylor2( @twoders, 0.0, 2.0, 0.5, 10);

the matrix A will contain the values of the arguments t_i 
and the approximate values of the solution w_i 
of Example 1 in Section 5.3 (page 268), 
and, if you type 

>> A

and hit "return", you will see values 
like in the corresponding column 
of Table 5.3 on page 270.