
function [x,t]=modEuler(vf,a,b,x0,N)
%solves dx/dt=vf(x,t) x(a)=x0 by explicit modified Euler  
%scheme xn=xp+h*vf(xp+0.5*hvf(xp,tp)) i.e. explic Runge scheme of 2nd order
%INPUT
%vf- function handle to vector field : function y=vf(x,t)... end
%a,b,- interval ends
%x0 - initial value in a column i.e. x(a)=x0
%N- nmber of discrete time steps -1 
%OUTPUT
%t - discrete time steps in vector t=[a, a+h, a+2*h,..,a+N*h]
%    h=(b-a)/N
%x - matrix N+1 x dim(x0) in (k+1)th row we have approximation of x(a+k*h)
%    computed by explicite Euler
%EXAMPLE:
%vf=@(x,t) x; N=10;a=0;b=1;x0=1;[x,t]=modEuler(vf,a,b,x0,N);plot(t,x,t,exp(t),";sol;")
%2D example:
%vfp=@(x,t) [0,1;-1,0]*x;sol=@(t)[cos(t);-sin(t)];
%N=40;a=0;b=8.3;x0=[1;0];[x,t]=modEuler(vfp,a,b,x0,N);plot(t,x,t,sol(t)');norm(x-sol(t)')

 h=(b-a)/N;%step
 t=linspace(a,b,N+1); %discrete times
 x=zeros(N+1,length(x0));
 x(1,:)=x0';xn=x0;
 for k=1:N,
  xp=xn;
  fp=vf(xp,t(k));
  xn=xp+h*vf(xp+0.5*h*fp,t(k)+0.5*h); %mod Euler
  %xn=xp+h*fp; %explicit Euler
  x(k+1,:)=xn';
 end%for k
end%function Taylor()
