%% physical parameters
a = 3; %u(x=0)=a
b= 3; %u(x=1) = b

%% set the parameters
N = 100; %number of interior nodes
h=1./(N+1);
x_array = 0:h:1; %define the physical domain, x = (0,1)
x_array = transpose(x_array(2:N+1 )); %remove the 0 and 1, which are the boundaries

%% construct the K matrix
K_mat = zeros(N, N);
for i=1:N
    for j=1:N
        if i==j
            K_mat(i,j) = 2./h^2;
        elseif abs(i-j)==1
            K_mat(i,j) = -1/h^2;
        end
    end
end
%% construct the right-hand side, f
f_vec =  cos(2.*x_array);

%% construct the additional vectors to take into account the boundary conditions
a_vec = zeros(N, 1);
a_vec(1) = a/h^2;
a_vec(N) = b/h^2;

%% solve for u
u_vec = K_mat\(f_vec  + a_vec);
u_analy =cos(2 * x_array)/4 - (cos(2) / 4 - b -1. /4 + a) * x_array -1./4 + a;


figure(1)
plot(x_array,u_analy)
hold on
plot(x_array,u_vec)
title("Analytical and Numerical Solution for -u''=cos(2x)")



%% Convergence Analysis
err_array = zeros(1,99);
for N=2:100
  %  N = 100; %number of interior nodes
    h=1./(N+1);
    x_array = 0:h:1; %define the physical domain, x = (0,1)
    x_array = transpose(x_array(2:N+1 )); %remove the 0 and 1, which are the boundaries

    %% construct the K matrix
    K_mat = zeros(N, N);
    for i=1:N
        for j=1:N
            if i==j
                K_mat(i,j) = 2./h^2;
            elseif abs(i-j)==1
                K_mat(i,j) = -1/h^2;
            end
        end
    end
    %% construct the right-hand side, f
    f_vec =  cos(2.*x_array);

    %% construct the additional vectors to take into account the boundary conditions
    a_vec = zeros(N, 1);
    a_vec(1) = a/h^2;
    a_vec(N) = b/h^2;

    %% solve for u
    f_vec(1) = f_vec(1)+a/h^2;
    f_vec(N) = f_vec(N)+b/h^2;
    u_vec = K_mat\(f_vec );

    u_analy = cos(2 * x_array)/4 - (cos(2) / 4 - b -1. /4 + a) * x_array -1./4 + a ;
    err_array(1,N-1) = sqrt(sum((u_vec-u_analy).^2)/N); %compute the L^2 error of the solution
end


figure(2)
fit_params = polyfit(log(2:1:100),log(err_array),1); %linear fit on log-log scale
slope = fit_params(1);
y_int = fit_params(2);


loglog((2:1:100),err_array ,'o')
hold on
loglog((2:1:100), (2:1:100).^slope * exp(y_int),'r')
title("L2 Norm Decay of the Solution for -u''=cos(2x)")
xlabel("N")
ylabel("L^2 Error")
txt = ['Slope = ' num2str(slope)  ];
text(10,0.0001,txt,'FontSize',14)