%-----------------------------------------------------%
% First order differential eq: modified Euler         %
%-----------------------------------------------------%
%INPUT: t1 - initial time
%       tN - final time
%       x1 - initial condition
%       N  - number of points on the grid
%       f  - derivative function (anonymous f(x,t))
%OUTPUT: xt_modE and t
 function [xt_modE, t] = modEuler_method(t1, x1, tN, N, f)

t=linspace(t1,tN,N);
dt=t(2)-t(1);

xt_modE(1)=x1;
for i=1:N-1
    xint(i)=xt_modE(i)+dt*f(xt_modE(i),t(i));
    xt_modE(i+1)=xt_modE(i)+dt*(f(xt_modE(i),t(i))+f(xint(i),t(i+1)))/2;
end

end