function CDteken
global np xy

set(gca,'xlim',[0 1],'ylim',[-0.1 1.1]);

clear x y xp xsort punt
hold on
for i=1:np,
   punt(i)=plot([xy(i,1) xy(i,1)] , [0  0] , 'kx' , ...
		'LineWidth',3, 'MarkerSize',12);
   eval(['set( punt(' , num2str(i) , '),''ButtonDownFcn'',' ...
	 ' ''pt=ginput(1);' ...
	 '   xy(', num2str(i), ',1)=pt(1);' ...
	 '   clf, CDteken'' )']);
end


% bepaal roosterpunten uit punten op scherm
xp(1)=0;
for i=1:np, 
  xp(i+1)=xy(i,1);
end
xp(np+2)=1;
xsort=sort(xp);
xp=xsort;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


N=np+1;

clear x xrek yexrek xrek0 yexrek0 b hcrs h Nc Nf 
clear M2L rl2L y2L M2S rl2S y2S upp low dia


%%%%%%%%%%%%%%%%%%%%%% grid generation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xrek(1)=xp(2);
h(1)=xrek(1);
for i=2:N,
  xrek(i)=xp(i+1);
  h(i)=xrek(i)-xrek(i-1);
end

u=1;
k=0.01;
rechts=0;

%%%%%%%%%%%%%%%%%%%%%%% determine exact solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% exacte oplossing op gerekt rooster 
%
pe=1/k;
% aktie om underflow van exp(-pe) te voorkomen

if k < 0.1,
hcrs= (1-10*k)/50;
hfin= 10*k/50;
  xex= [0: hcrs : 1-10*k   1-10*k+hfin : hfin : 1]; 
else
  xex=[0:0.01:1];
end

if pe > 100,
   empe=0;
else
   empe=exp(-pe);
end
yex=1-(1-exp(pe*(xex-1)))/(1-empe); % exact op fijn rooster

yexrek=1-(1-exp(pe*(xrek-1)))/(1-empe); % exact op rekenrooster
%
b = zeros(1,N-1);

xrek0=[0 xrek];
yexrek0=[0 yexrek];
yexrek=yexrek';
 
%%%%%%%%%%%%%%%%%%%%  2e orde methoden  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Tweede orde Lagrange discretisatie
%
for i=1:N-1
   upp(i) = (h(i)* u - 2*k) / h(i+1);
   low(i) = (-h(i+1) * u - 2*k) / h(i);
   dia(i) = -upp(i)-low(i);
   rl2L(i) = b(i) * (h(i)+h(i+1));   
end
%
% matrix
M2L=diag(upp(1:N-2),1) + diag(dia(1:N-1),0) + diag(low(2:N-1),-1);
%
% boundary conditions
rl2L(N-1)=rl2L(N-1)-upp(N-1);
%
% solution
y2L = M2L\rl2L';
y2L(N) = 1;
if N==2, y2L=y2L'; end % Matlab kan niet zien of y2L een rij- dan wel
                       % kolomvector was


%
% 2nd order symmetry-preserving method **********
%
for i=1:N-1
   upp(i) = u - 2*k/h(i+1);
   low(i) = -u - 2*k/h(i);
   dia(i) = -upp(i)-low(i);
   rl2S(i) = b(i) * (h(i)+h(i+1));
end
%
% matrix
M2S=diag(upp(1:N-2),1) + diag(dia(1:N-1),0) + diag(low(2:N-1),-1);
% 
% boundary conditions
rl2S(N-1)=rl2S(N-1)-upp(N-1);
%
% solution
y2S = M2S\rl2S';
y2S(N) = 1;
if N==2, y2S=y2S'; end


% alles plotten

pSiser=exist('pS');
if pSiser==1,
   set(pS(1),'visible','off');
   set(pS(2),'visible','off');
   set(pS(3),'visible','off');
end
pS=plot(xex,yex,'g',xrek0,[0 y2L'],'r', xrek0,[0 y2S'],'b');
set(pS(1),'hittest','off');
set(pS(2),'hittest','off');
set(pS(3),'hittest','off');
set(pS(1),'Linewidth',3);
set(pS(2),'Linewidth',3);
set(pS(3),'Linewidth',3);
set(gca,'XLim',[0 1],'Ylim',[-0.1 1.1]);
[abc def]=legend('exact','Lagrange','symmetric',2);
set(abc,'Fontsize',16);
xd=get(def(3),'xdata');
yd=get(def(3),'ydata');
set(def(3),'xdata',[xd-0.1 xd+0.1],'ydata',[yd yd]);
set(def(3),'linestyle','-','linewidth',3,'color','g','marker','none');
xd=get(def(5),'xdata');
yd=get(def(5),'ydata');
set(def(5),'xdata',[xd-0.1 xd+0.1],'ydata',[yd yd]);
set(def(5),'linestyle','-','linewidth',3,'color','r','marker','none');
xd=get(def(7),'xdata');
yd=get(def(7),'ydata');
set(def(7),'xdata',[xd-0.1 xd+0.1],'ydata',[yd yd]);
set(def(7),'linestyle','-','linewidth',3,'color','b','marker','none');