function [x,y]= odeeuler(fvec, xspan, y0, N)
	% solves system of ordinary differential equation with Euler's method
	% x -  column vector of x positions
	% y -  solution array values of y, each row corresponds to  particular row of x.
	% fvec - handle to a function f(x,y) returning forces column vector
	% xspan vector with initial and final x coordinates i.e. [x0, xf]
	% y0 - initial conditions
	% N - number of points in the x column (N>=2)

	x0=xspan(1);
	xf=xspan(2);



	h=(xf-x0)/(N-1); % step
	x=linspace(x0,xf,N); % values of x

	odeorder=length(y0);
	y=zeros(N,odeorder); % initialization

	x(1)=x0; y(1,:)=y0;  % initial conditons
	for i=2:N % number of steps is less by 1 then number of points since we know x0,y0
		xprev=x(i-1);
		yprev=y(i-1,:);
		% matlab somehow always send column vector for 'y' thus 
		% transposing yprev
		f=fvec(xprev, yprev.'); 
		% to make this method compatible with matlab
		% we expect f as a column vector 
		% thus we need to transpose 
		f=f.'; % note the dot in .' this avoid complex conjugate transpose
		ynext=yprev+f*h;  % vector of new values of y: y(x+h)=y(x)+f*h
		y(i,:)=ynext;
	end
end




