echo on 
#####################################################################
#
#
#		Proste przyklady na rozwiazywanie 
#                  rownan rozniczkowych zwyczajnych
#
#####################################################################

#Zdefuniujmy pole wektorowe (skalarne)
###############################################################33
function y=pw(x,t)
y=-x;
end

#Rozwiazemy rownanie x'=-x =pw(x) z x(0)=1 na odcinku t=0,..,1
# na siatce z h=0.1
#i porownamy z rozwaizaniem dokladnym (w pktach siatki)
t=0:0.1:1;
sol=exp(-t); #rozwiazanie dokladne w pktach siatki 
x=lsode("pw",1,t);
norm(x-sol',"inf")

#teraz wezmy 2 wymiarowe rownanie x''=-x x(0)=0 x'(0)=1
#rozwiazanie sin(x)
#Definijemy pole wektorowe x'=y, y'=-x
function y=pw(x,t)
y(1)=x(2);
y(2)=-x(1);
end
t=0:0.1:7;
x=lsode("pw",[0 1],t);
#Pierwsza kolumna to rozwiazanie druga jego pochodna
#blad rozwiazania
 norm(x(:,1)- sin(t)',"inf")
#i jego pochodnej
norm(x(:,2)- cos(t)',"inf")
#I zrobmy wykres 
plot(t,sin(t),";rozw dokladne;",t,x(:,1),"g+;rozw przyblizone;");
#Pokrywa sie!
#A dla chetnych wykres trajektorii x versus x' - powinien wyjsc okreg
#plot(x(:,1),x(:,2));
#plot3(x(:,1),x(:,2),t);

