

function [y,t]=midpoint(vf,t0,x0,x1,n,t)
#Function solves approx. dx/dt=vf(x,t);x(t0)=x0
#using midpoint scheme i.e. .computes x(tkt0+k*h) k=0,..,n
#Output
#y=[y1 ... y(n+1)] -> y(k+1) approximate sol(t0+k*h)
#t=[t0, t0+h...t0+n*h]
#Input
#vf  - function vf(x,t)- vector field - with input x and t and ouput value of vf(x,t)
#t0, x0 values from initial value
#x1 approximation of x(t0+h) for h=(t-t0)/n
#n- number of steps 
#t - time value
h=(t-t0)/n;
t=zeros(1,n+1);t(1)=t0;t(2)=t0+h;
y=zeros(length(x0),n+1);
y(:,1)=x0;y(:,2)=x1;

for k=3:(n+1),
 y(:,k)=y(:,k-2)+2*h*vf(y(:,k-1),t(k-1));
 t(k)=t(k-1)+h;
endfor
endfunction

