-----------------------------------------------------
Example of Matlab Code 	
-----------------------------------------------------

Solve the IVP for ODE: y' = f(y,x),  y(x0) = y0

A) Euler method (EULER.m): 

function y = euler(func, x0, y0, xn, h)

% The inputs func, x0 and y0 define an IVP y' = func(x,y), y(x0) = y0.
% The file "EULER.m" provides an approximation to y(xn) using Euler method 
% with step size h.

% initialize
f = inline(func);
x = x0;
y = y0;
n = (xn-x0)/h;

% iterate
for i = 1:n
	k1 = f(x, y);
	y = y + h*k1;
	x = x + h;
end


B) Runge-Kutta method (RK.m): 

function y = rk(func, x0, y0, xn, h)

% The inputs func, x0 and y0 define an IVP y' = func(x,y), y(x0) = y0.
% The file "RK.m" provides an approximation to y(x0:h:xn) using classical 
% fourth-order Runge-Kutta method with step size h.

% initialize
f = inline(func);
x(1) = x0;
y(1) = y0;
n = (xn-x0)/h;

% iterate
for i = 1:n
	k1 = f(x(i), y(i));
	k2 = f(x(i) + h/2, y(i) + h/2*k1);
	k3 = f(x(i) + h/2, y(i) + h/2*k2);
	k4 = f(x(i) + h, y(i) + h*k3);
	y(i+1) = y(i) + h/6*(k1 + 2*k2 + 2*k3 + k4);
	x(i+1) = x(i) + h;
end

C) Predictor-Corrector method (PC.m): 

function y = pc(func, x0, y0, xn, h)

% The inputs func, x0 and y0 define an IVP y' = func(x,y), y(x0) = y0.
% The file "PC.m" provides an approximation to y(x0:h:xn) using predictor-corrector
% method with fourth-order Adams-Bashforth (predictor) and fourth-order
% Adams-Moulton (corrector) using step size h.

% initialize
f = inline(func);
x(1) = x0;
y(1) = y0;
n = (xn-x0)/h;

% generate starting estimates using Runge-Kutta
for i = 1:3
	k1 = f(x(i), y(i));
	k2 = f(x(i) + h/2, y(i) + h/2*k1);
	k3 = f(x(i) + h/2, y(i) + h/2*k2);
	k4 = f(x(i) + h, y(i) + h*k3);
	y(i+1) = y(i) + h/6*(k1 + 2*k2 + 2*k3 + k4);
	x(i+1) = x(i) + h;
end

% iterate
for i = 4:n
	% Adams-Bashforth -- *predict*
	y(i+1) = y(i) + h/24*(55*f(x(i), y(i)) - 59*f(x(i-1),y(i-1)) ...
		+ 37*f(x(i-2),y(i-2)) - 9*f(x(i-3),y(i-3)));
	x(i+1) = x(i) + h;
	% Adams-Moulton -- *correct*
	y(i+1) = y(i) + h/24*(9*f(x(i+1),y(i+1)) + 19*f(x(i),y(i)) ...
		- 5*f(x(i-1),y(i-1)) + f(x(i-2),y(i-2)));
end


COMPARE the results (compare.m): 

function compare(varargin)
% COMPARE.M
%
% Display a table containing data sets for different methods. The inputs
% should be in pairs data, description. It is assumed that x = 0:.1:1, and
% this is displayed as the first column in the resulting table.
%
% Author: Luis Carvalho (carvalho@dam.brown.edu)
% Last update: 03/05/04

nsets = nargin/2;

hbase = 0.1;
for i = 1:nsets
	data{i} = varargin{2*i-1};
	desc{i} = varargin{2*i};
	h(i) = 1/(length(data{i})-1);
	step(i) = hbase/h(i);
end

x = 0:hbase:1;
for i = 1:nsets
	y{i} = data{i}(1:step(i):length(data{i}));
end

%header
header = '\tx';
for i = 1:nsets
	header = strcat(header, '\t', desc{i}, '\t');
end
disp(sprintf(header));
%body
for i = 1:11
	body = '\t%0.1f';
	for j = 1:nsets
		body = strcat(body, '\t%1.8f');
		output(j) = y{j}(i);
	end
	disp(sprintf(body, x(i), output));
end