%% Diferenční metoda pro lineární okrajový problém
% obecná rovnice: p(x)y'' + q(x)y' + r(x)y=f(x)
% Příklad: y''=y/x^2-y'/x, y(1)=1, y(5)=0.2, přesné řešení y=1/x
% Definice funkcí:

p=@(x) ones(size(x));
q=@(x) 1./x;
r=@(x) -1./(x.^2);
f=@(x) zeros(size(x));

% Další parametry - interval, okrajové podmínky a počet subintervalů:

I=[1 5];
BV=[1 0.2];
n=5;
[x1,y1]=solveBP1(p,q,r,f,I,BV,n);
plot(x1,y1,x1,1./x1,'--')
title('Okrajovy problem y''''=y/x^2-y''/x, y(1)=1, y(5)=0.2');
legend('numericke reseni, n=5','presne reseni')
fprintf('Max. chyba, 5 subintervalu: %f\n',max(abs(y1-1./x1)));

%

n=10;
[x1,y1]=solveBP1(p,q,r,f,I,BV,n);
plot(x1,y1,x1,1./x1,'--')
title('Okrajovy problem y''''=y/x^2-y''/x, y(1)=1, y(5)=0.2');
legend('numericke reseni, n=10','presne reseni')
fprintf('Max. chyba, 10 subintervalu: %f\n',max(abs(y1-1./x1)));

%

n=100;
[x1,y1]=solveBP1(p,q,r,f,I,BV,n);
plot(x1,y1,x1,1./x1,'--')
title('Okrajovy problem y''''=y/x^2-y''/x, y(1)=1, y(5)=0.2');
legend('numericke reseni, n=100','presne reseni')
fprintf('Max. chyba, 100 subintervalu: %f\n',max(abs(y1-1./x1)));

%
%% Diferenční metoda pro lineární okrajový problém
% speciální rovnice: -(p(x)y')' + q(x)y = f(x)
% Příklad: (y'/x)'=4*x*y, y(0.1)=exp(-0.01), y(2)=exp(-4), přesné řeěení: y(x)=exp(-x^2)
% Definice funkcí:

p=@(x) 1./x;
q=@(x) 4*x;
f=@(x) zeros(size(x));

% Další parametry - interval, okrajové podmínky a počet subintervalů:

I=[0.1 2];
BV=[exp(-0.01) exp(-4)];
n=5;
[x1,y1]=solveBP2(p,q,f,I,BV,n);
plot(x1,y1,x1,exp(-x1.^2),'--')
title('Okrajovy problem (y''/x)''=4*x*y, y(0.1)=exp(-0.01), y(2)=exp(-4)');
legend('numericke reseni, n=5','presne reseni')
fprintf('Max. chyba, 5 subintervalu: %f\n',max(abs(y1-exp(-x1.^2))));

%

n=10;
[x1,y1]=solveBP2(p,q,f,I,BV,n);
plot(x1,y1,x1,exp(-x1.^2),'--')
title('Okrajovy problem (y''/x)''=4*x*y, y(0.1)=exp(-0.01), y(2)=exp(-4)');
legend('numericke reseni, n=10','presne reseni')
fprintf('Max. chyba, 10 subintervalu: %f\n',max(abs(y1-exp(-x1.^2))));

%

n=100;
[x1,y1]=solveBP2(p,q,f,I,BV,n);
plot(x1,y1,x1,exp(-x1.^2),'--')
title('Okrajovy problem (y''/x)''=4*x*y, y(0.1)=exp(-0.01), y(2)=exp(-4)');
legend('numericke reseni, n=100','presne reseni')
fprintf('Max. chyba, 100 subintervalu: %f\n',max(abs(y1-exp(-x1.^2))));

%

function [x,y] = solveBP1(p,q,r,f,I,BV,n)
% function [x,y] = solveBP1(p,q,r,f,I,BV,n)
% Solving the linear boundary problem by Finite Differences Method
% The equation: p(x)y'' + q(x)y' + r(x)y=f(x)
% p,q,r,f   - annonymous functions
% I=[a,b]   - interval
% BV        - boundary values, BV(1)=y(a), BV(2)=y(b)
% n         - number of subintervals, equally-spaced knots are generated
%
% x - knots
% y - approximate solution in the knots
%
% Finite Differences:   y''(x) ~ [y(x-h)-2y(x)+y(x+h)]/h^2
%                        y'(x) ~ [y(x+h)-y(x-h)]/2h
%

ya=BV(1);yb=BV(2);  % boundary values
a=I(1);b=I(2);
h=(b-a)/n;
x=linspace(a,b,n+1);
x=x(:);     % column
xx=x(2:n);  % inner points
P=p(xx);
Q=q(xx);
R=r(xx);
F=f(xx);

d=-2*P/h^2+R;                       % main diagonal of the 3-diagonal system
l=P(2:end)/h^2-Q(2:end)/(2*h);      % subdiagonal
u=P(1:end-1)/h^2+Q(1:end-1)/(2*h);  % superdiagonal
b=F;                                % the right side of the system
b(1)=b(1)-P(1)*ya/h^2+Q(1)*ya/(2*h);            % y(a) is moved to the right side
b(end)=b(end)-P(end)*yb/h^2-Q(end)*yb/(2*h);    % y(b) is moved to the right side

% solving of the 3-diagonal system follows

yy=zeros(n-1,1);  % solution in inner points

% Gaussian elimination
for k=1:n-2
    q=l(k)/d(k);
    % l(k)=0;   % zeroing of the subdiagonal, not necessary
    d(k+1)=d(k+1)-q*u(k);
    b(k+1)=b(k+1)-q*b(k);
end

% back substitution
yy(n-1)=b(n-1)/d(n-1);
for k=n-2:-1:1
    yy(k)=(b(k)-yy(k+1)*u(k))/d(k);
end

y=[ya;yy;yb];

end     % That's all, Folks!

function [x,y] = solveBP2(p,q,f,I,BV,n)
% function [x,y] = solveBP2(p,q,f,I,BV,n)
% Solving the linear boundary problem by Finite Differences Method
% The equation: -(p(x)y')' + q(x)y = f(x)
% !!! Be careful of the sign '-' on the previous line. !!!
%
% p,q,f     - annonymous functions
% I=[a,b]   - interval
% BV        - boundary values, BV(1)=y(a), BV(2)=y(b)
% n         - number of subintervals, equally-spaced knots are generated
%
% x - knots
% y - approximate solution in the knots
%
% Finite Differences: 
% -(p(x)y'(x))' ~ -[p(x-h/2)y(x-h)-(p(x-h/2)+p(x+h/2))y(x)+p(x+h/2)y(x+h)]/h^2
%

ya=BV(1);yb=BV(2);  % boundary values
a=I(1);b=I(2);
h=(b-a)/n;
x=linspace(a,b,n+1);
x=x(:);                     % column
xx=x(2:n);                  % inner points
x1=(x(1:end-1)+x(2:end))/2; % midpoints of the subintervals
P=p(x1);
Q=q(xx);
F=f(xx);

d=(P(1:end-1)+P(2:end))/h^2+Q;  % main diagonal of the 3-diagonal system
l=-P(2:end-1)/h^2;              % subdiagonal
%u=-P(2:end-1)/h^2;             % superdiagonal
u=l;                            % superdiagonal is the same as l
b=F;                            % the right side of the system
b(1)=b(1)+P(1)*ya/h^2;          % y(a) is moved to the right side
b(end)=b(end)+P(end)*yb/h^2;    % y(b) is moved to the right side


% solving of the 3-diagonal system follows

yy=zeros(n-1,1);  % solution in inner points

% Gaussian elimination
for k=1:n-2
    q=l(k)/d(k);
    % l(k)=0;   % zeroing of the subdiagonal, not necessary
    d(k+1)=d(k+1)-q*u(k);
    b(k+1)=b(k+1)-q*b(k);
end

% back substitution
yy(n-1)=b(n-1)/d(n-1);
for k=n-2:-1:1
    yy(k)=(b(k)-yy(k+1)*u(k))/d(k);
end

y=[ya;yy;yb];

end     % That's all, Folks!