%% Metoda střelby
% Rovnice y''=f(x,y,y'), y(a)=y_a,y(b)=y_b
% Řešíme počáteční problém s y'(a)=\alpha, hledáme takovou hodnotu \alpha, že platí druhá okrajová podmínka, tj. řešíme rovnici F(\alpha)=0,\ \mathrm{pro}\ F(\alpha)=y(b,\alpha)-y_a.
% Rovnici převádíme na systém rovnic
% y_1'=y_2, y_2'=f(x,y_1,y_2), y_1(a)=y_a, y_2(a)=\alpha
% Příklad: y''=-\frac{1}{4xy}, y(1)=1,y(4)=2, přesné řešení y=\sqrt x.
% Rovnice:

f=@(x,y)[y(2);-1./(4*x.*y(1))];

% Další parametry a počáteční iterace pro \alpha.

format short
a=1;b=4;
I=[a,b];
y10=1;
y1n=2;
al0=(y1n-y10)/(b-a)
y0=[y10;al0];

% Pro numerické řešení použijeme program ode45.

[xx,yy]=ode45(f,I,y0);

% Hodnota F(\alpha_0):

F0=yy(end,1)-y1n

%
% Metoda sečen - provedeme korekci \alpha a určíme hodnotu funkce F:

al1=al0-F0/(b-a)
y0=[y10;al1];
[xx,yy]=ode45(f,I,y0);
F1=yy(end,1)-y1n

%
% Další iterace \alpha metodou sečen:

al2=al1-F1*(al1-al0)/(F1-F0)
y0=[y10;al2];
[xx,yy]=ode45(f,I,y0);
F2=yy(end,1)-y1n
al3=al2-F2*(al2-al1)/(F2-F1)
y0=[y10;al3];
[xx,yy]=ode45(f,I,y0);
F3=yy(end,1)-y1n
al4=al3-F3*(al3-al2)/(F3-F2)
y0=[y10;al4];
[xx,yy]=ode45(f,I,y0);
F4=yy(end,1)-y1n
format long
al5=al4-F4*(al4-al3)/(F4-F3)
y0=[y10;al5];
[xx,yy]=ode45(f,I,y0);
F5=yy(end,1)-y1n
plot(xx,yy(:,1),xx,sqrt(xx),'--')

%
% Příklad 2: lin. rovnice y''=\frac{1}{4x^2}y-\frac{y'}{x}, y(1)=1,y(4)=2, přesné řešení y=\sqrt x.
% Rovnice:

f=@(x,y)[y(2);y(1)./(4*x.^2)-y(2)./x];

% Další parametry a počáteční iterace pro \alpha.

format short
a=1;b=4;
I=[a,b];
y10=1;
y1n=2;
al0=(y1n-y10)/(b-a)
y0=[y10;al0];

% Pro numerické řešení použijeme program ode45.

[xx,yy]=ode45(f,I,y0);

% Hodnota F(\alpha_0):

F0=yy(end,1)-y1n

%
% Provedeme korekci \alpha a určíme hodnotu funkce F:

al1=al0-F0/(b-a)
y0=[y10;al1];
[xx,yy]=ode45(f,I,y0);
F1=yy(end,1)-y1n

%
% Další iterace \alpha metodou sečen:

al2=al1-F1*(al1-al0)/(F1-F0)
y0=[y10;al2];
[xx,yy]=ode45(f,I,y0);
F2=yy(end,1)-y1n

% Protože se jednalo o lin. rovnici, druhá iterace už je prakticky přesná.