clear;

for i=1:5

n=10*(2^(i-1))-1; h=1/(n+1);

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

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

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

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

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

%plot((0:n+1)'*h,sol,'b')
%hold on;
%plot((0:n+1)'*h,exact,'r')
%hold off;

hoo(i)=h;
err1(i)=max(abs(sol-exact));
est(i)=(10^4/96)*h^2;

clear('sol','sol1','exact')

end
