
function [x,t]=eulerz17(pw,x0,t0,T,n=100)
%prosta implementacja niejawnego/(zamknietego) schematu Eulera
%implicit Euler
%INPUT
%pw pole wektorowe /vector field
%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 implicit Euler i.e. x(k+1)=x(k)+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
%[x,t]=eulerz17(pw=@(t,x) x,x0=1,t0=0,T=1,n=100);plot(t,x,";zE;",t,exp(t),";exp(t);");
%2D example
%graphics_toolkit('gnuplot')
%[x,t]=eulerz17(pw=@(t,x) [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),";zE - 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]=eulerz17(pw=@(t,x) [x(2);-sin(x(1))],x0=[0;1],t0=0,T=7,n=100);
%plot(x(:,1),x(:,2),";zE - 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,
%explicit Euler  x(:,k+1)=x(:,k)+h*pw(t(k),x(:,k));
 %nonlinear function to be solved in order to get next step approx: 
 %x(k+1)-h*f(t(k+1),x(k+1)-x(k)=0
 nf=@(y) y-h*pw(t(k+1),y)-x(:,k);
 %initial starting point x(k)-the approx. sol. at the previous time step t(k) 
 x(:,k+1)=fsolve(nf,x(:,k));
 end
 x=x'; %to stay consistent with a standard octave ODE solver lsode()
end
