############################################
# Predyktor korektor zamkn/otwarty Euler
#  nieliniowy uklad rozwiazuje funkcja fsolve
############################################


#definiujemy zmienne globalne potrzebne do funkcji
#ktorej zerem jest kolejne przyblizenie x(tn) w zamknietym schemacie Eulera
global xn h tn



function y=f(x,t)
 #pole wektorowe
 y(1)=x(2);
 y(2)=-sin(x(1));
 y=y';
end

function y=G(x)
 #funkcja ktorej zerem jest x_{n+1} dla zadanych h n i x_n
 global xn h tn 

 y=x-h*f(x,tn+h) - xn;
end

 #Definiujemy ilosc krokow na start
 N=100;

 h=2*pi/N;

 x=zeros(2,N+1);
 t=zeros(1,N+1);
 x(:,1)=[1 ;0];
 t(1)=0;
 tn=t(1);

#petla calkujaca zamk Euler
 for k=1:N,
  xn=x(:,k);
  y=xn+h*f(xn,tn);
  x(:,k+1)=fsolve("G",y);
  tn=tn+h;
  t(k+1)=tn;
 endfor  

#A teraz porownamy z lsode (oczywiscie lsode ma nieporownywalnie lepszy rzad...
X=lsode("f",x(:,1),t);
h*norm(X-x',1) #dzielimy przez N aby uniezaleznic sie od N w normie
#a na koniec wykres trajektrii x vs y
opis=[";zamkn E h=" num2str(h) ";"];
plot(x(1,:),x(2,:),opis,X(:,1),X(:,2),";lsode;");
#jest kiepsko ale krok zamknietego Eulera dosc spory jak na schemat rzedu 1!

