%% set the parameters
a = 0;
b = 0;
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 = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY

%% solve for u
u_vec = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY

%% Analytical solution
u_analy = ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY


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


%% 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 =  ;%  ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY

    %% solve for u
    u_vec =  ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY

    %% Analytical solution
    u_analy =  ;% ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIFY
    
    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)