function [yvec,tvec] = midpoint(fun,a,b,h,y0)
% Task : midpoint to solve an ODE


t = a;
y = y0;
tvec = [a];
yvec = [y0];

while t+h < b
	yprime1 = feval(fun, [t; y]);
	yprime = feval(fun, [t+h/2,y+h/2*yprime1]);
	t = t + h;
	y = y + h*yprime;
	tvec = [tvec t];
	yvec = [yvec y];
end


return;