function [x,t]=euler(fun,x0,t0,tk,n)
%[x,t]=erluer(fun,x0,t0,tk,n)
%finds an approximate solution the differential equation
%x'(t)=fun(t,x(t))
%x(t0)=x0
%on the interval [t0,tk] using the Eulers method 
%n is the number of steps
d=length(x0);%the dimension of the system (number of equations)
t=linspace(t0,tk,n);%independent variable
h=t(2)-t(1);%step size
x=zeros(d,n);%space for the solution
x(:,1)=x0;%initial value 
for i=1:n-1
	k1=h*feval(fun,t(i),x(:,i));
	x(:,i+1)=x(:,i)+k1;
end