#zad 1 rano
disp("kolokwium ON gr 830 - maj 2008");
function [x  cnd]=SVAXF(N)
#funkcja rozwiazujaca uklad z Ax=eN z A NxN taka ze A(i,i)=4 
# A(i,i+1)=A(i+1,i)=1 i=1,..,N-1 
#A(N,1:N-1)=A(1:N-1,N)=1 
#i wektorem prawej strony eN N-tym wersorem


I=J=V=zeros(N+2*(N-1)+2*(N-2));
#Gl przek Naddiag Poddiag ostatni rzad ost kol
I=[1:N 1:(N-1) 2:N N*ones(1,N-2) 1:(N-2)];
J=[1:N 2:N 1:(N-1) 1:(N-2)       N*ones(1,N-2) ];
V=[4*ones(1,N) ones(1,4*N -6)];
A=sparse(I,J,V,N,N);
eN=zeros(N,1);
eN(N)=1;
x=A\eN;

#sprawdzmy czy macierz dobrze rozwiazuje....
if(norm(A*x-eN,"inf")>1e-6) 
 printf("Blad residualny podejrzanie duzy!\n"); 
endif
#bez tego if zawsze funkcja liczylaby cond co jest drogie nawet 
#jak uzytkownik chce tylko rozwiazanie
 if(nargout>1)
  cnd=cond(A);
 endif
# if(N<10) 
#  full(A) # to tylko do testu czy macierz OK
# endif
endfunction

#testujemy czy działa 
[x  cnd]=SVAXF(10)

# a teraz porownajmy czasy wywolan
tic;SVAXF(1000);toc
#i z obl cond
tic;[x cnd]=SVAXF(1000);toc
#ostrzezenie stad ze macierz jest symetryczna osobliwa 
#ale ma ujemne wart wlasne - a octave probuje ja rozwiazac met Choleskiego - ktora #wtedy nie dziala

#zad 2 - wykres trajektorii na [0,1] dla rownania [x;y]'=[-1 308;-308 -1]*x;
#z [x;y](0)=[1;1]
global a 
a=308
function y=pw(x)
global a
y=[-1 a;-a -1]*x;
end


t=linspace(0,1);
X=lsode("pw",[1;1],t);
plot(X(:,1),X(:,2),";x vs y;");pause(2);
#oj a malo punktow ewidentnie
t=linspace(0,1,1e3);
X=lsode("pw",[1;1],t);
plot(X(:,1),X(:,2),";x vs y;");pause(2);
#jeszcze za malo 
t=linspace(0,1,1e4);
X=lsode("pw",[1;1],t);pause(2);
plot(X(:,1),X(:,2),";x vs y;")
