clear all
close all
% -u''(x)+x*u(x)=(1+x)*cos(x)
% u(0) = 1 (u'(0) = 0)
% u(pi/2) = 0
a = 0;
b = pi/2;
m = 100;
h = (b-a)/(m-1);
x = linspace(a,b,m)';
q = x;
g = (x+1).*cos(x);
D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
% D2(1,1:2) = [-2,2]/h^2; % Neumann
A = -D2+diag(q);
A(1,1:2) = [1,0]/h^2; % Dirichlet
%M = 2^60; % penalizzazione
%A(1,1) = M; % penalizzazione
A(m,m-1:m) = [0,1]/h^2;
%A(m,m) = M; % penalizzazione
g(1) = 1/h^2; % Dirichlet
%g(1) = M; % penalizzazione
g(m) = 0;
u = A\g;
% check if A is spd
R=chol(A);
norm(R'*R-A,'fro') % Frobenius norm is easy to compute
plot(x,u,'*')
xlabel('x')
ylabel('u')
% convergenza
mrange = 10:10:100;
counter = 0;
for m = mrange+1
  counter = counter+1;    
  h = (b-a)/(m-1);
  x = linspace(a,b,m)';
  q = x;
  g = (x+1).*cos(x);
  D2 = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
%  D2(1,1:2) = [-2,2]/h^2; % Neumann
  A = -D2+diag(q);
  A(1,1:2) = [1,0]/h^2; % Dirichlet
  A(m,m-1:m) = [0,1]/h^2;
  g(1) = 1/h^2; % Dirichlet
  g(m) = 0;
  u = A\g;
  errore(counter) = norm(u-cos(x),inf);
end
figure
loglog(mrange,errore,'*',mrange,mrange.^(-2))
legend('errore','ordine 2')
xlabel('m')
ylabel('e_m')
