% finite difference method

h = 0.25;
x_vec = 0 : h : 1;
n_x = length(x_vec);

f = @(x) 2;

%c_vec = zeros(n_x-2,1);
f_vec = zeros(n_x-2,1);

for i = 1:n_x-2
    
    xim1 = x_vec(i);
    xi = x_vec(i+1);
    xip1 = x_vec(i+2);
    q = @(x) ((x-xim1)/h)*f(x);
    p = @(x) ((xip1 -x)/h)*f(x);
    
    f_vec(i) = integral(q,xim1,xi) + integral(p,xi,xip1);
    
    
end


L_h = zeros(n_x-2,n_x-2);
for i = 1:n_x -2
    for j = 1:n_x -2
        
        if i == j 
            L_h(i,j) = 2/h;
        elseif i == j +1
            L_h(i,j) = -1/h;
        elseif j == i +1
            L_h(i,j) = -1/h;
        end
        
        
    end
end




c_vec = L_h\(f_vec);


u_sol = @(x) -x^2 + x;

%u_num = @(x)  phi_coe(1,1)*(subplus(x_vec()) - 2*subplus(x-1) + subplus(x-2));

hold on

fplot(u_sol, [0,1],'Color','black');


plot(x_vec,[0,c_vec',0]);
legend('u_{sol}', 'u_{num}');


hold off