function func = fd1ord ( imeth, iout ),
% FD1Ord
%
% This m-file is part of the material of the lectures on
%   Numerical Solution of Partial Differential Equations
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  May 10, 2006
%
% Plot output data for semi-infinite string
%
% Parameters:
%   imeth  (input)  : control parameter method
%                       =1 ... Explicit Euler
%                       =2 ... Lax-Friedrichs
%                       =3 ... Upwind
%                       =4 ... Lax-Wendroff
%   iout   (input)  : control parameter output
%                       =0 ... output of numerical solution
%                       =1 ... output of error for selected h, \tau
%                       =2 ... output of error for selected h and \tau/h=const
%
% Examples:
%   for imeth=1:4, for iout=0:2,
%      fd1ord ( imeth, iout );  
%      eval(sprintf('print -dpng ../png/fd1ord%1i%1i.png',imeth,iout));
%   end; end;

% -> control data plot
graphwidth   = 14.0;
graphheight  = 14.0;
fac	     =  0.9;
shift	     =  0.5;
border       =  2.0;

ActFontSize  = 14;
ActLineWidth =  2;

%
% -> Bestimme Groesse des Bildschirms und Skalierungsfaktor fuer die Ausgabe
root = 0;
set ( root, 'Units', 'centimeters' );
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);
set ( root, 'Units', 'pixels' );
screenscale = 1.0;

%
% -> Definiere den Bildschirmausschnitt
if length ( get ( root, 'Children' ) ) == 0,
  h = figure;
else
  h = gcf;
  clf;
end;
stdunits = get ( h, 'Units' );
set ( h, 'Units',    'centimeters' );
set ( h, 'Position', ...
  [ left+width-border-screenscale*graphwidth height-border-screenscale*graphheight screenscale*graphwidth screenscale*graphheight ] );
set ( h, 'Units', stdunits );

%
% -> Skalierung der Ausgabe in PostScript
%
set ( gcf, 'PaperUnits',    'centimeters' );
set ( gcf, 'PaperPosition', [ 0 0 graphwidth graphheight ] );

% control parameters output
nxx = 41;
ntt = 11;
xxmin = -2.0;
xxmax =  2.0;
ttmin =  0.0;
ttmax =  1.0;
xx = linspace ( xxmin, xxmax, nxx );
tt = linspace ( ttmin, ttmax, ntt );

umin =  0.0;
umax =  1.2;

errmax = [ 0.02 0.1 0.02 5e-4 ];
tstr = [ '(Euler)         '
         '(Lax-Friedrichs)'
         '(Upwind)        '
         '(Lax-Wendroff)  ' ];

% define the problem and the space and time grid
c = 1.5;

nx = 5*([ 501 1001 501 1001 ]-1)+1;
nt =    [ 201  201 401  401 ];
xmin = -24.0;
xmax =  26.0;
tmin =  ttmin;
tmax =  ttmax;

% plot 4 diagrams for iout>0 and only one diagram for iout=0
if iout~=0,
   minplt = 1;
else
   minplt = 4;
end;

% -> prepare output
clf,

for iplot=minplt:4,

   % -> stepsizes
   if iout~=2,
      nxact = nx(iplot);
      ntact = nt(iplot);
   else
      nxact = 2^(iplot-1)*(nx(1)-1)+1;
      ntact = 2^(iplot-1)*(nt(1)-1)+1;
   end;
     
   h   = (xmax-xmin)/(nxact-1);
   tau = (tmax-tmin)/(ntact-1);
   lam = tau/h;

   e = ones ( nxact, 1 );
   if imeth==1,
      a = speye(nxact) - lam*c/2 * spdiags([-e e],[-1 1],nxact,nxact);
   elseif imeth==2,
      a = spdiags([e e],[-1 1],nxact,nxact)/2 - lam*c/2 * spdiags([-e e],[-1 1],nxact,nxact);
   elseif imeth==3,
      if c<0,
         a = speye(nxact) - lam * c * spdiags([-e e], [0 1],nxact,nxact);
      elseif c>=0,
         a = speye(nxact) - lam * c * spdiags([-e e], [-1 0],nxact,nxact);
      end;
   elseif imeth==4,
      a = speye(nxact) - lam*c/2 * spdiags([-e e],[-1 1],nxact,nxact) + ...
          lam^2*c^2 / 2 * spdiags([ e -2*e e],-1:1,nxact,nxact);
   end;

   % -> initialize workspace
   x = linspace ( xmin, xmax, nxact );
   t = linspace ( tmin, tmax, ntact );
   u = zeros ( nxact, ntact );
   u(:,1) = exp(-x.^2)';
   
   % -> numerical solution
   for it=1:(ntact-1),
      u(:,it+1) = a * u(:,it);
   end;
  
   % -> compute error
   delx = round ( (nxact-1) / (nxx-1) * (xxmax-xxmin) / (xmax - xmin) );
   delt = (ntact-1) / (ntt-1);
   ixl  = round ( (nxact-1) * (xxmin-xmin) / (xmax-xmin) + 1 );
   ixr  = round ( (nxact-1) * (xxmax-xmin) / (xmax-xmin) + 1 );

   [ XX, TT ] = meshgrid ( xx, tt );
   err = abs ( u(ixl:delx:ixr,1:delt:ntact) - exp(-(XX-c*TT).^2)' );

   if iout==0,
      % -> plot the solution
      mesh ( tt, xx, u(ixl:delx:ixr,1:delt:ntact) );
      set ( gca, 'YDir', 'reverse' );
      axis ( [ ttmin, ttmax, xxmin, xxmax, umin, umax ] );
      xlabel ( 't' );
      ylabel ( 'x' );
      zlabel ( 'u(x,t)' );
      title ( ['u_t + 1.5 u_x = 0,    u(x,0) = exp(-x^2)     ', tstr(imeth,:) ] );
      view ( -37.5-25, 20 );
   else
      % -> plot error
      subplot ( 2, 2, iplot );
      mesh ( tt, xx, err );
      set ( gca, 'YDir', 'reverse' );
      set ( gca, 'ZLim', [ 0 errmax(imeth) ] );
      xlabel ( 't' );
      ylabel ( 'x' );
      tstr1 = sprintf('h = 1/%3i,', round(1/h) );
      tstr2 = sprintf(' = 1/%3i', round(1/tau) );
      title ( [ tstr1, ' \tau', tstr2 ] );
      view ( -37.5-25, 35 );
   end;

end;

%
% -> Korrektur der Groesse und Lage der Koordinatensysteme, Aenderung der Schriftgroesse
allax = get ( gcf, 'Children' )';
for ax = allax,
  rect    = get ( ax, 'Position' );
  rect(2) = rect(2) + shift * (1-fac) * rect(4);
  rect(4) = fac * rect(4);
  set ( ax, 'Position', rect );
  set ( ax, 'FontSize', ActFontSize );
  set ( get ( ax, 'XLabel' ), 'FontSize', ActFontSize );
  set ( get ( ax, 'YLabel' ), 'FontSize', ActFontSize );
  set ( get ( ax, 'ZLabel' ), 'FontSize', ActFontSize );
  set ( get ( ax, 'Title' ), 'FontSize', ActFontSize );
end;
