MATLAB lessons N.3

A simple implementation of Euler scheme: File Euler.m

:open:
 1%
 2% This routine approximate the solution of the IVP (Initial Value Problem)
 3%
 4%  /  y' = f(t,y)
 5%  |
 6%  |  y(a) = y[0]
 7%   \
 8%
 9%  PARAMETERS:
10%   FUN = reference to user defined function
11%   t0  = Starting integration time
12%   y0  = Initial data of the IVP
13%   h   = time step
14%   N   = number of step to be performed
15%
16%  OUTPUT
17%   T   = vector containing the nodes where 
18%         the approximate solution is computed
19%         t = t0, t0+h, t0+3*h,..... t0+N*h
20%   Y   = vector containing the approximate solution at nodes T, i.e
21%         Y(1) ~ y(t0), 
22%         Y(2) ~ y(t0+h),
23%         Y(3) ~ y(t0+2*h),
24%         .....
25%         Y(N+1) ~ y(t0+N*h)
26%
27function [T,Y] = Euler( FUN, t0, y0, h, N )
28  % pre allcation of the output vector T and Y
29  T = zeros(N+1,1) ;
30  Y = zeros(N+1,1) ;
31  
32  % set T(1) and Y(1) with the initial condition
33  T(1) = t0 ;
34  Y(1) = y0 ;
35  
36  % loop performing the Euler Steps
37  for k=1:N
38    T(k+1) = T(k) + h ;
39    Y(k+1) = Y(k) + h * feval( FUN, T(k), Y(k) ) ;
40    % feval( FUN, T(k), Y(k) ) == FUN( T(k), Y(k) )
41  end
42
43end

A simple ODE to test the code: File problem1.m

:open:
1%
2% function for the IVP
3%
4%  y' = -t*y
5%
6function dy = problem1( t, y )
7  dy = -t*y ;
8end