odeeuler.m

function [t,y] = odeeuler(func,a,b,y0,h);
% solution routine for Euler's method for solving
%
%      dy/dt = func(t,y)
%
% USAGE:  [t,y] = odeeuler(func,a,b,y0,h)
%
% input  func    = name of external function to evaluate the RHS
%                  of the ODE
%        a, b    = limits of integration
%        y0      = initial condition
%        h       = stepsize
%
% output [t, y]  = solution vectors

t = [a];
y = [y0];
i = 1;

while t(i) < b

   i = i+1;

   t(i) = a + (i-1)*h;

   y(i) = y(i-1) + h*feval(func,t(i-1),y(i-1));

end