
function [x,t]=taylor19(f,ft,fx,x0,t0,T,n=100)
%prosta implementacja otwartego schematu Taylora rz 2
%INPUT
%f,ft,fx pole wektorowe /vector field f(x,t) and its derivatives f_t f_x
%x0 warunek poczatkowy/initial value
%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 Tyalorem tj x(k+1)=x(k)+h*f(t(k),x(k)) +0.5*h*h...;
%     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 Taylor scheme
%     x(k) approximates x(t(k)) x - the solution of IVP
%t - discrete times t=[t0,t0+h,..T]
%Przyklad/Example
%n=20;T=1;[x,t]=taylor19(f=@(x,t) x,ft=@(x,t) 0,fx=@(x,t) 1,x0=1,t0=0,T,n);plot(t,x,";Taylor;",t,exp(t),";exp(t);");
%n=20;T=10;[x,t]=taylor19(f=@(x,t) -x,ft=@(x,t) 0,fx=@(x,t) -1,x0=1,t0=0,T,n);plot(t,x,";Taylor;",t,exp(-t),";exp(-t);");

 x=zeros(length(x0),n+1);
 t=linspace(t0,T,n+1)';
 h=T/n;
 x(:,1)=x0;
 for k=1:n,
  fk=f(x(:,k),t(k));
  x(:,k+1) =x(:,k)+h*fk;
  x(:,k+1)+=0.5*h*h*(ft(x(:,k),t(k))+fx(x(:,k),t(k))*fk);
 end
 x=x'; %to stay consistent with a standard octave ODE solver lsode()
end
