clear
xdel(winsid());

x0 = 0;
y0 = 1;
xe = 2;
ye = 2;

Nx = 5;
x = linspace(x0,xe,Nx)';
dx = x(2)-x(1);

function wyn = fa(x)
    wyn = -sqrt(5)*x/6 + sqrt(1 + x.^2)/3 + (2+x.^2)/3;
endfunction

function wyn = f1(x)  // elements of the main diagonal
     wyn = 4*(1+x.^2)+2*dx^2;
endfunction

function wyn = f2(x) // elements of the upper diagonal
     wyn = 2*(1+x.^2)+x*dx;
endfunction

function wyn = f3(x) // elements of the lower diagonal
     wyn = 2*(1+x.^2)-x*dx;
endfunction

function wyn = f4(x)   //elements of the right han d side vector
     wyn = 2*x.^2*dx^2;
endfunction

A = diag(-f1(x(2:$-1))) + diag(f2(x(2:$-2)), 1) + diag(f3(x(3:$-1)),-1);

b = f4(x(2:$-1));
b(1) = b(1) - f3(x(2))*y0;
b($) = b($) - f2(x($-1))*ye;

y = A\b;
y=[y0;y;ye];

ya = fa(x);
scf();

plot2d(x,y,style=5)
plot2d(x,ya,style=2)
xtitle(msprintf("Numerical and analytical aolutions, %d division",Nx));
legend(['numerical','analytical'],2,%f);

scf();
plot2d(x,y-ya,style=5)
xtitle(msprintf("Absolute errors, %g divisions",Nx));