
function [x,t]=midpoint17(f,x0,x1,t0,T,n=100)
%prosta implementacja otwartego schemat umidpoint -2krokowego
%INPUT
%f pole wektorowe /vector field f(x,t)
%x0 warunek poczatkowy/initial value
%x1 - starting value aproximating x(t0+h) computed by a one-step scheme in practise
%t0 start czasu/inital time
%T - koniec odcinka calkowania
%n - ilosc pktow czasowych tj krok czasowy h=T/n
%OUTPUT
% x - wektor dlugosci n+1 z przyblizeniem rozwiazania dx/dt=pw(t,x) x(t0)=x0
%     obliczonymi midpointem tj x(k+2)=x(k)+2*h*pw(t(k+1),x(k+1);
%     x(k) przybliza x(t(k)) x - rozwiazanie dokladne
%     vector of length n+1 approximating solution of IVP: dx/dt=pw(t,x) x(t0)=x0
%     COMPUTED BY midpoint i.e. x(k+2)=x(k)+2*h*pw(t(k+1),x(k+1);
%     x(k) approximates x(t(k)) x - the solution of IVP
%t - discrete times t=[t0,t0+h,..T]
%Przyklad/Example
%n=100;t0=0;T=1;[x,t]=midpoint17(pw=@(x,t) -x,x0=1,x1=exp(-(T-t0)/n),t0,T,n);plot(t,x,";midp;",t,s=exp(-t),";exp(-t);");pause(2);plot(x-s,";error;");
 x=zeros(length(x0),n+1);
 t=linspace(t0,T,n+1)';
 h=T/n;
 x(:,1)=x0;
 x(:,2)=x1;
 for k=1:n-1,
  fkp1=f(x(:,k+1),t(k+1));
  x(:,k+2)=x(:,k)+2*h*fkp1;
 end
 x=x'; %to stay consistent with a standard octave ODE solver lsode()
end
