

1==1;
#solver octave'a lsode

#najprostsze wywolanie
#dla y'=y y(0)=1

function y=f(x,t)
 y=x;
endfunction

x0=1;
t0=0;
t=[t0,1];
X=lsode(@f,x0,t)

#x(k+1) przyblizone rozwiazanie dla czasu z wektora t=[t0,t1,..,tN] 
#tzn dla tk u nas t[t0=0,t1=1] czyli X(2) przybliza x(1)=exp(1)
X(2)-exp(1) 

#wykres dla y'=ay y(t0)=x0 (jak podac a do pola wektorowego?
global a
function y=f(x,t)
 global a
 y=a*x;
endfunction
x0=1;
a=1;T=1;
t=linspace(0,T);
X=lsode(@f,x0,t);
plot(t,X,";lsode;",t,exp(a*t),";exp(at);");
pause(3);
#moze dla wiekszego a i dluzszego czasu
x0=1;
a=10;T=2;
t=linspace(0,T);
X=lsode(@f,x0,t);
S=exp(a*t)';
plot(t,X,";lsode;",t,S,";exp(at);");
pause(2);
#a blad?
bb=X-S;
bw=bb./S;
max([bb,bw])
#bezwgledny wyglada strasznie ale wzgledny nie taki zly....

# rownanie wahadla zlinearyzowane

function y=f(x,t)
#r wahadla zlinearyzowane y''=-y
 y=zeros(2,1);
 y(1)=x(2);
 y(2)=-x(1);
endfunction

#dla 4 pkt czasowych
X=lsode(@f,[0;1],pi*[0,0.25,0.5,1])
#wykres
x0=[0;1];;
a=1;T=10;
t=linspace(0,T);
X=lsode(@f,x0,t);
plot(t,X(:,1),";wahadlo y;",t,X(:,2),";y';");
pause(1);
plot(X(:,1),X(:,2),";trajek y vs y';");
pause(1);


#zaleznosc do parametru
#y'=y*(1-a*y)
#y(0)=y0
#rozwiazanie zalezy od s i a tzn y(t,s,a)
#narysuj wykres zaleznosci od s dla t=1 i a=0 i 1
#dla s na [-0.5,3] tzn F0(s)=y(1,s,0)
function y=f(x,t)
global a
 y=x*(1-a*x);
end
s=linspace(-.5,3);
global a
a=0;
F0=zeros(length(s),1);
for k=1:length(s),
 X=lsode(@f,s(k),[0,1]);
 F0(k)=X(2);
endfor
plot(s,F0,";s->y(1,s) y=y, y(0)=s;");
pause(1)
#prosta nie zaskakuje znamy rozwiazanie
#y(t,s,0)=s*exp(t) czyli F0(s)=s*exp(1)
a=1;
F1=zeros(length(s),1);
for k=1:length(s),
 X=lsode(@f,s(k),[0,1]);
 F1(k)=X(2);
endfor
plot(s,F1,";s->y(1,s) y=y(1-y), y(0)=s;");
pause(1);
#sprobujcie Panstwo policzyc rozwiazanie dla s=-1....

#a teraz wykres wzgledem parametru a na [-0.5,1] z t=1 s=1
aa=linspace(-0.5,1);
s=1;
T=1;
Fa=zeros(length(aa),1);
for k=1:length(aa),
 a=aa(k);
 X=lsode(@f,s,[0,T]);
 Fa(k)=X(2);
endfor
plot(aa,Fa,";a->y(1,a) y=y(1-ay), y(0)=1;");
pause(1)

#prosze sprobowac policzyc wartosc y(1,s=1,a=-1)
#dlaczego sie nie udaje i  lsode() drukuje ostrzezenia?
#oczywiscie jako zadanie pozostwiam policzenie analityczne
#i wtedy wszystko staje sie jasne.... mozna uzyc maxima'e 




