% FEM 1D problem pg. 96

close all 
clear all

c0 = 299792458;
eps0 = 8.8541878e-12;
mu0 = pi*4e-7;


alpha = 1/mu0;
beta = -omega^2*eps0;

a = 2e-2;       % Glass thickness (=2cm)
omega = 3*10^9;
E0 = 1;         % amplitude of incident field E^i_z
er = 2.5;       % relative permittivity 
sigma = 0.02;   % conductivity

% Boundary conditions
fa = 1;
fb = 0;

% Subdivide domain into 7 elements (8 nodes)    
i = 1:8; % node index
x = i-3; % coordinate for each node
a = x(1);
b = x(end);
f = zeros(size(i));
f(1) = fa;
f(end) = fb;



basis_funcs = zeros(8,8);
for idx = 1:8
    basis_funcs(idx, idx) = 1;
end

%w =  % weights. Galerkins method dictates that they equal the basis functions. We exclude bounds since the sought field is known there (Dirichlet).

basis_funcs_ddx = diff(basis_funcs, 1, 2);
phi = basis_funcs;
dphi = basis_funcs_ddx;

ksi = @(i,j) sum(alpha * dphi(i,j) .* dphi(i,j) + beta * phi(i,j) .* phi(i,j));

A = spdiags([1, 1, 1],-1:1, 6, 6);
for i = 2:6
    for j = i-1:i+1
        if j < 2 || j > 6
        else
            A(i,j) = ksi(i,j);
        end
    end
end

b = zeros(6, 1);
b(1) = -fa*ksi(2,1);
b(end) = fb;

z = A\b;
z = [fa;z;fb]

plot(x, z)

% figure 
% hold on
% grid on
% ylim([-2, 3])
% for idx = i
%     plot(x, basis_funcs(i,:))
% end