expeuler.m

ode routines in matlab

function [x, y] = ExpEuler( x0, x_end, y0, h, rhs_fctn )
%
%  Purpose:
%  Solve the initial value problem
%               
%                      y'(x) = f(x,y(x))
%                      y(x0) = y0
%
%  on the interval [ x0, x_end ] using the explicit Euler method.
%  The explicit Euler method is given by
%
%    y(:,i+1) = y(:,i) + h * f(x(i),y(:,i))
%
%
%  Usage:
%
%       [x, y] = ExpEuler( x0, x_end, y0, h, rhs_fctn )
%
%  Parameters:
%
%  Input:
%
%  x0          initial x
%
%  x_end       final x
%
%  y0          initial values
%
%  h           step size
%
%  rhs_fctn    function evaluating f(x,y) for given x and y
%              rhs_fctn(x, y, '') returns  f(x,y)
%
%
%
%  Output:
%  
%  x        vector of x values x(i) = x0 + (i-1)*h
%
%  y        y(:,i) contains approximation of solution at x(i)
%

% find number of time steps needed to reach final time
mx = ceil((x_end-x0)/ h);


% allocate workspace
n = size(y0(:),1);
x = (x0:h:x0+mx*h);
y = zeros(n,mx+1);

% Set initial values
y(:,1) = y0(:);

for i = 1:mx
    y(:,i+1) = y(:,i) + h * feval(rhs_fctn,x(i),y(:,i),'');
end