% To Do: f = x^2

clear all
close all

N = 100;
h = 1/N;
a = 2;
b = -1;
c = -1;

x = linspace(0,1,N+2);

A = diag(a*ones(1,N)) + diag(b*ones(1,N-1),1) + diag(c*ones(1,N-1),-1);

b = ones(N,1);
xi = (1/h)*A\b;

sol(1) = 0;
sol(2:N+1) = xi;
sol(N+2) = 0;

plot(x,sol,'x-','linewidth',2)



return
clear all
close all

%
%Domain: [0,1]; nx=100;
%dx=1/(nx-1);f(x)=delta(1/2)


%
nx = 10;
dx = 1/nx;


% source term
%Boundary conditions:
%u(0)=u(1)=0
s=(1:nx)*0;s(nx/2)=1.;
u1=0; s(1) =0;
unx=0; s(nx)=0;

% Boundary conditions:
% u(0)=0.15
% u(1)=0.05
s=(1:nx)*0;s(nx/2)=1.;
u1=0.15; s(2) =u1/dx;
unx=0.05; s(nx-1)=unx/dx;
% 

% assemble matrix Aij
A=zeros(nx);
for i=2:nx-1
    for j=2:nx-1
        if i==j
            A(i,j)=2/dx;
        elseif j==i+1
            A(i,j)=-1/dx;
        elseif j==i-1
            A(i,j)=-1/dx;
        else
            A(i,j)=0;
        end
    end
end
fem(2:nx-1)=inv(A(2:nx-1,2:nx-1))*s(2:nx-1)';
fem(1)=u1;
fem(nx)=unx;


plot(fem)

