clear;
n=10; h=1/(n+1);

xdisc=((0:(n+1))')*h;
exact=exsol(xdisc);

%% S1

mdiag1=diag(ones(n+1,1));
ldiag1=diag(ones(n,1),-1);
udiag1=diag(ones(n,1),1);

amat1=(-ldiag1+2*mdiag1-udiag1);
amat1(n+1,n+1)=1;

xdisc1=((1:(n+1))')*h;
b1=h^2*rhsfu(xdisc1); b1(n+1)=h;

sol1=amat1\b1;
sol1=[0; sol1];

%% S2

mdiag2=diag(ones(n+2,1)); 
ldiag2=diag(ones(n+1,1),-1);
udiag2=diag(ones(n+1,1),1);

amat2=-ldiag2+2*mdiag2-udiag2;
amat2(n+2,n+2)=1;
amat2(n+2,n+1)=0;
amat2(n+2,n)=-1;

xdisc2=((1:(n+2))')*h;
b2=h^2*rhsfu(xdisc2); b2(n+2)=2*h;

sol2=amat2\b2;
sol2=[0; sol2(1:n+1)];

plot((0:n+1)'*h,sol1,'b')
hold on;
plot((0:n+1)'*h,sol2,'g')
plot((0:n+1)'*h,exact,'r')
hold off;

h
err1=max(abs(sol1-exact))
err2=max(abs(sol2-exact))

