% the program takes the intial data y_0 at the initial time x_0 and 
% computes up to time x_final using N timesteps of size dx.  It creates
% two vectors, x and y, which you can then plot against each other.
% The ODE you're trying to solve is dy/dx = f(x,y) where the vector
% field is contained in the subroutine f.m
%
% [x,y] = euler_step(y_0,x_0,x_final,N,'f')

function [x,y] = euler_step(y_0,x_0,x_final,N,f)

% for dy/dx = f(y(x),x)
% y(x+dx) ~ y(x) + dx y'(x) = y(x) + dx f(y(x),x)

% step-size:
dx = (x_final-x_0)/N;

% put initial data into vector
x(1) = x_0;
y(1) = y_0;
for ii=1:N
  F1 = feval('f',x(ii),y(ii));
  y(ii+1) = y(ii) + dx*F1;
  x(ii+1) = x_0 + ii*dx;
end