x = 0: 0.1: 4;

deff('y = f(x)','y = (16 - x^2)^(1/2)');
y = f(x);

//draw
xgrid();
plot2d(x, y);

k = splin(x, y);
t = 0: 1: 4;
spl = interp(t, x, y, k);
plot2d(t, spl, -10);
plot2d(f(1),1,-10);

t(6) = t(5);
spl(6) = spl(5);
spl(5) = 1.0;
t(5) = 3.87;

//table of coordinates
disp("x   -   y");
for i = 1:1:length(t)
    mprintf('%2.2f - %2.2f\n', t(i), spl(i));
end

//u(x, y)
function u = Calc_U(x, y)
    u = 0.5 * abs(x) + abs(y);
endfunction
for i = 1:1:length(t)
    u(i) = Calc_U(t(i), spl(i));
end

disp("U(x, y) = 0.5 * |x| + |y|");
for i = 1:1:(length(t))
    mprintf('U(%5.2f, %5.2f) = %5.2f\n', t(i), spl(i), u(i));
end

//Matrix A
A = [ 4 -2  0 -1  0  0  0  0  0  0  0  0  0  0  0;
     -1  4 -1  0 -1  0  0  0  0  0  0  0  0  0  0;
      0 -1  4  0  0 -1  0  0  0  0  0  0  0  0  0;
     -1  0  0  4 -2  0  0 -1  0  0  0  0  0  0  0;
      0 -1  0 -1  4 -1  0  0 -1  0  0  0  0  0  0;
      0  0 -1  0 -1  4 -1  0  0 -1  0  0  0  0  0;
      0  0  0  0  0 -1  4  0  0  0 -1  0  0  0  0;
      0  0  0 -1  0  0  0  4 -2  0  0 -1  0  0  0;
      0  0  0  0 -1  0  0 -1  4 -1  0  0 -1  0  0;
      0  0  0  0  0 -1  0  0 -1  4 -1  0  0 -1  0;
      0  0  0  0  0  0 -1  0  0 -1  4  0  0  0 -1;
      0  0  0  0  0  0  0  0 -4  0  0  4  0  0  0;
      0  0  0  0  0  0  0  0 -2  0  0 -1  4 -1  0;
      0  0  0  0  0  0  0  0  0 -2  0  0 -1  4 -1;
      0  0  0  0  0  0  0  0  0  0 -1  0  0 -1  4;];
B = [4; 4.37; 8.92; 0; 0; 0; 8.3; 0; 0; 0; 2.94; 0; 0; 0; 0;];

C=rref([A B]);
[n,m]=size(C);
x=C(:,m)

disp(x);