function integrate(mode)

global panel_discr panel_timestep panel_tfin panel_points panel_diff 
global panel_grid panel_stretch plotfig

clear x xrek yexrek xrek0 yexrek0 b hcrs h Nc Nf L 
clear N grid method k delt tfin 
clear M2L uppL diaL lowL rl2L y2L y2l y2loud V2L D2L ewL
clear M2S uppS diaS lowS rl2S y2S y2s y2soud V2S D2S ewL

N=str2num(get(panel_points,'String'));
grid=get(panel_grid, 'Value');      % 1 = uniform, 2 = stretched 3 = abrupt
method=get(panel_discr,'Value');    % 1 = method A, 2 = method B, 3 = both
L=str2num(get(panel_stretch,'String'));

k=str2num(get(panel_diff,'String'));
delt=str2num(get(panel_timestep,'String'));
tfin=str2num(get(panel_tfin,'String'));

%N=4;
%grid=3;
%method=2;
%L=3;
%k=0.01;
%delt=0.01;
%tfin=5.0;


% convective velocity
u=1;

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

Nc=round(0.5*N);
Nf=N-Nc;

if grid==2 & L*k < 0.5,
% stretched grid
   fac=(L*k/(1-L*k))^(1/Nc);
%   disp(['stretch factor= ',num2str(fac)])
   mesh=(1-fac)/(1-fac^N);
   xrek(1)=mesh;
   for i=2:N
     mesh=fac*mesh;
     xrek(i)=xrek(i-1)+mesh;
   end
elseif grid==3 & L*k < 0.5,
% abrupt grid
   hcrs=(1.-L*k)/Nc;           % coarse grid cells left
   hfin=L*k/Nf;                 % fine grid cells right
   xrek(1:Nc)=hcrs*(1:1:Nc);                 % grid points left
   xrek(Nc+1:N)=xrek(Nc) + hfin*(1:1:Nf);   % grid points right
else
% uniform grid
   xrek(1:N)=(1:1:N)/N;
end

%xrek

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

%h


%%%%%%%%%%%%%%%%%%%%%%% determine exact solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% exacte solution on fine grid with 100 grid points
%
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

pe=1/k;
% action to prevent underflow in exp(-pe)
if pe > 100,
   empe=0;
else
   empe=exp(-pe);
end
yex=1-(1-exp(pe*(xex-1)))/(1-empe); % exact on fine grid

yexrek=1-(1-exp(pe*(xrek-1)))/(1-empe); % exact on selected grid
%

xrek0=[0 xrek];
yexrek0=[0 yexrek];
yexrek=yexrek';
 

%%%%%%%%%%%%%%%%%%%%  2nd order methods  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 2nd order Lagrange discretisation
%
for i=1:N-1
   uppL(i) = (h(i)* u - 2*k) / h(i+1);
   lowL(i) = (-h(i+1) * u - 2*k) / h(i);
   diaL(i) = -uppL(i)-lowL(i);
   rl2L(i) = 0;   
end
%
% matrix
M2L=diag(uppL(1:N-2),1) + diag(diaL(1:N-1),0) + diag(lowL(2:N-1),-1);
%
% boundary conditions
rl2L(N-1)=rl2L(N-1)- uppL(N-1);
%
% solution
y2L = M2L\rl2L';
y2L(N) = 1;
if N==2, y2L=y2L'; end % Matlab cannot see whether y2L is a row or 
                       % a column vector 


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


if N==2, y2S=y2S'; end



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

if (mode==1), 

% plot steady solution

figure(3)
clf

if (method == 2)  ,    % show only Lagrange
  pS= plot(xex,yex,'g', xrek0,[0 y2L'],'r-*');
  set(gca,'XLim',[0 1],'Ylim',[-0.1 1.1]);
  set(pS(1),'Linewidth',3);
  set(pS(2),'Linewidth',3);
  title('Steady solution','fontsize',16)

  % legend maken
  [abc def]=legend('exact','Lagrange',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');

else

  pS= plot(xex,yex,'g', xrek0,[0 y2L'],'r-*', xrek0,[0 y2S'],'b-');
  set(gca,'XLim',[0 1],'Ylim',[-0.1 1.1]);
  set(pS(1),'Linewidth',3);
  set(pS(2),'Linewidth',3);
  set(pS(3),'Linewidth',3);
  title('Steady solution','fontsize',16)

  % legend maken
  [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');

end % plot modus

end % mode=1


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


if mode==2,

% timeintegration

t=0;
y2l=xrek';
y2s=xrek';

figure(2)
clf
set(gca,'XLim',[0 1],'Ylim',[-1.0 2.0]);

while t < tfin,

   t=t+delt;

% Lagrange
   y2loud=y2l;
   i=1;
   y2l(i) = y2loud(i) - delt * (diaL(i)*y2loud(i) + uppL(i)*y2loud(i+1));
   for i=2:N-1,
      y2l(i) = y2loud(i) - delt*(lowL(i)*y2loud(i-1) + diaL(i)*y2loud(i) ...
                                                   + uppL(i)*y2loud(i+1));
   end
   y2l(N)=1;

% symmetry-preserving
   y2soud=y2s;
   i=1;
   y2s(i) = y2soud(i) - delt*(diaS(i)*y2soud(i) + uppS(i)*y2soud(i+1));
   for i=2:N-1,
      y2s(i) = y2soud(i) - delt*(lowS(i)*y2soud(i-1) + diaS(i)*y2soud(i) ...
                                                   + uppS(i)*y2soud(i+1));
   end
   y2s(N)=1;

   if (method == 2),                    % show only Lagrange 
      ps1=plot(xrek0, [0 y2l'],'r-*');
      set(gca,'XLim',[0 1],'Ylim',[-1.0 2.0]);
      set(ps1(1),'Linewidth',3);
      tijdstring=sprintf('%6.3f', t);
      title(['time  ', tijdstring],'Fontsize',20);
   else
      ps=plot(xrek0, [0 y2l'],'r-*', xrek0, [0 y2s'],'b');
      set(gca,'XLim',[0 1],'Ylim',[-1.0 2.0]);
      set(ps(1),'Linewidth',3);
      set(ps(2),'Linewidth',3);
      tijdstring=sprintf('%6.3f', t);
      title(['time  ', tijdstring],'Fontsize',20);
   end % plot modus
   drawnow


end % time step

end % mode=2


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

if (mode == 3)   % show eigenvalues

   [V2L,D2L]=eig(M2L);   
   [V2S,D2S]=eig(M2S);

for i=1:N-1
   ewL(i)=-D2L(i,i);
   ewS(i)=-D2S(i,i);
end

%close(figure(4))
figure(4)
clf

if (method==2),  % show only Lagrange 

  plot(0,0,'.', real(ewL), imag(ewL),'r*', 'markersize',12,'linewidth',2)
  xgrens = get(gca,'XLim');
  ygrens = get(gca,'YLim');
  hold on
  plot([0  0], [ygrens(1) ygrens(2)],'k','linewidth',2)
  plot([xgrens(1) xgrens(2)],[0 0], 'k','linewidth',1)
  title('Eigenvalues','fontsize',16)

  % make legend 
  [abc2 def2]=legend('Lagrange',2);
  set(abc2,'Fontsize',16);
  set(def2(3),'color','r','marker','*','markersize',12,'linewidth',2);

else

  plot(real(ewL), imag(ewL),'r*', real(ewS), imag(ewS),'bo', ...
      'markersize',12,'linewidth',2)
  xgrens = get(gca,'XLim');
  ygrens = get(gca,'YLim');
  hold on
  plot([0  0], [ygrens(1) ygrens(2)],'k','linewidth',2)
  plot([xgrens(1) xgrens(2)],[0 0], 'k','linewidth',1)
  title('Eigenvalues','fontsize',16)

  % make legend 
  [abc2 def2]=legend('Lagrange','symmetric',2);
  set(abc2,'Fontsize',16);
  set(def2(3),'color','r','marker','*','markersize',12,'linewidth',2);
  set(def2(5),'color','b','marker','o','markersize',12,'linewidth',2);

end % plot modus

end % mode=3