
function [x,t]=eulero17(f,x0,t0,T,n=100)
%prosta implementacja otwartego schematu Eulera
%INPUT
%f pole wektorowe /vector field f(x,t)
%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 otwartym Eulerem tj x(k+1)=x(k)+h*pw(t(k),x(k);
%     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 EXPLICIT Euler i.e. x(k+1)=x(k)+h*pw(t(k),x(k);
%     x(k) approximates x(t(k)) x - the solution of IVP
%t - discrete times t=[t0,t0+h,..T]
%Przyklad/Example
%[x,t]=eulero17(f=@(x,t) x,x0=1,t0=0,T=1,n=100);plot(t,x,";oE;",t,exp(t),";exp(t);");
%2D example
%graphics_toolkit('gnuplot')
%[x,t]=eulero17(pw=@(x,t) [x(2);-x(1)],x0=[0;1],t0=0,T=7,n=100);
%plot(t,x(:,1),";oE;",t,sin(t),";sin(t);");pause(2);plot(x(:,1),x(:,2),";oE - trajectory;",sin(t),cos(t),";sol;");
%x(n+1,1)-sin(t(n+1))%as we know the solution...
%really  nonlinear  2D ODE - pendulum ODE
%[x,t]=eulero17(pw=@(x,t) [x(2);-sin(x(1))],x0=[0;1],t0=0,T=7,n=100);
%plot(x(:,1),x(:,2),";oE - trajectory;");
%we do not know the exact solution...
 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;
 end
 x=x'; %to stay consistent with a standard octave ODE solver lsode()
end
