
% How to use diag
D=[1 2 3 4 5];
diag(D)

D=[1 2 3 4];
diag(D,1)

D=[1 2 3 4];
diag(D,-1)
%%
% Now a first example for finite difference approximations
% Let's solve f''=1; f(0)=f(1)=1; explicit solution: f(x)=1-x/2+x^2/2
N=100;
L=diag(-2*ones(1,N+1)) + diag(ones(1,N),1) + diag(ones(1,N),-1);
L(1,1) = 1;
L(1,2) = 0;
L(N+1,N+1) = 1;
L(N+1,N) = 0;

cond(L)

dx = 1/N;
b=ones(N+1,1)*dx^2;
A=1; B=1;
b(1) = A;
b(N+1) = B;

% the solution F approximates [f(0), f(1/N), ... f(1)]
F=L\b;

% now calculate the true F
X=linspace(0,1,N+1);
F_true=1-X/2+X.^2/2;
F_true=F_true';


figure;
plot(X,F,'*','MarkerSize',10)

hold on
plot(X,F_true,'r','LineWidth',4)

% maximum error

max(abs(F_true-F))

%%
% let's play with the right hand side

b(11) = b(11) + 5*dx^2;
F=L\b;
hold on
plot(X,F);

b(11) = b(11) + 50*dx^2;
F=L\b;
hold on
plot(X,F);





