clear();

function y = u_exact(x)
    y=x^(-1).*exp(x);
endfunction

function y=g_0(x)
    y =exp(2)/1.5;
  // y=exp(2)/2.*x^2-exp(2).*x;
endfunction

function y=g_k(k,x)
    y = -1 - cos(%pi*k*x);
  //  y = (x-1)^(k+1).*(2-x)^2;
endfunction

function y =r(x)
   // y=x.*(x+4);
   y=x;
endfunction

function y=p(x)
   // y=-(2*x+4);
   y=2;
endfunction

function y=q(x)
 //   y=ones(1,length(x))*(-x);
    y=-x;
endfunction

function y = Lfik(k,x)
    y = r(x).*cos(%pi*k*x)*(%pi*k)^2 + p(x).*sin(%pi*k*x)*%pi*k + q(x).*(1 - cos(%pi*k*x));
    //y = r(x).*(k*(k+1).*(x-1)^(k+1).*(2-x)^2 - (k+1).*(x-1)^k*2.*(2-x) -(k+1).*(x-1)^k*2.*(2-x) + (x-1)^(k+1).*2) + p(x).*((k+1).*(x-1)^k.*(2-x)^2 - (x-1)^(k+1).*2.*(2-x)) + q(x).*g_k(k,x);
endfunction

/*function y = gauss(A,b)
    N = length(b);
    for i = 1:N // Зануляемый столбец
        for j = 1:N
            A(i,1:M) = A(i,1:M)/A(i,j)
    end
endfunction
*/

h=0.01;
x = 1:h:2;
N = length(x);
//x_colloc = 0.005 + modulo(rand(1,N,'uniform'),( 1.995 - 0.005)); //  равномерное распределение на (0.005,1.995) - вектор из N эл-в
x_colloc = x(2:N-2);
M = length(x_colloc); // Кол-во узлов коллокации

A = zeros(M,M);
for k=1:M
    A(1:M,k) = Lfik(k,x_colloc); // Вставляю сразу столбцы коэффициентов при c_k
end


    b=-exp(2).*x_colloc'/1.5;
b=-3.89.*x_colloc';
//b=-exp(2)/2.*x_colloc'^2-exp(2).*x_colloc'.*x_colloc';
//b = -exp(2)/1.5 *ones(M,1); // Правая часть

c = linsolve(A,b);
y = g_0(x);
//y=3.89;
for k = 1:M
     y = y + c(k)*g_k(k,x);
end
//y=y-1.51;
u = u_exact(x);
xset('thickness',2);
plot(x, y,'b');
plot(x, u,'r');
legend ('y', 'u',2);
title('Collocation');