function [t,y]=sec1_4_burgersFDMFixedMesh(jmax)
%
% example driver for Burgers' equation (with exact solution)
% in section 1.4. it calls pdefdm() and lin_interp_err() and
% uses a fixed, uniform mesh.
%

%
% Copyright (C) 2010 Weizhang Huang and Robert D. Russell
% all rights reserved.
%
% This program is provided "as is", without warranty of any kind.
% Permission is hereby granted, free of charge, to use this program
% for personal, research, and education purposes. Distribution or use
% of this program for any commercial purpose is permissible
% only by direct arrangement with the copyright owner.
%

global epsilon;

  epsilon = 1.0e-4;

cpu0=clock;

   % job = 1 for solution and error at t = T
   % job = 3 for computing the global error.

   job = 1;

   %jmax = 41;
   npde = 1;
   nn = npde*jmax;

   x=zeros(jmax,1);
   u=zeros(npde,jmax);

% define initial solution

   a = 0;
   b = 1;
   x=linspace(a,b,jmax)';
   for j=1:jmax, u(1,j)=uexact(0,x(j)); end
   y=reshape(u,nn,1);
   yt = zeros(nn,1);
   
   reltol = 1e-6;
   abstol = 1e-4;
   
%   plot(x,u,'-o');
%   return;

% call the pde solver
   
   if job==1 % for solution and error at t = T
   
      tspan = [0 0.2 0.4 0.6 0.8 1.0];
      [t,y]=pdefdm(npde,jmax,tspan,y,yt,@PDE_F,@PDE_G,@BC_L,@BC_R,a,b,reltol,abstol);
      fprintf('\n cpu time used = %e \n', etime(clock,cpu0));
      % output the solution
      N=size(y,1)
	  t(N)
	  figure
	  labels={};
      marks='+osdv^<>';
      lines1='--:-';
      lines2='-  .';
      colors='mkbgrc';
      for n=1:N
         u=reshape(y(n,1:nn),npde,jmax);
         hold on
         labels{n} = ['t = ' num2str(t(n)) ];
         style=['-' marks(1+rem(n,8))];
         plot(x,u(1,:)',style);%,'LineWidth',2);
      end
      hold off
      legend(labels{:});
      xlabel('x');
      ylabel('u');
      axis([0 1 -0.2 1.2]);
	  box on;
	  
	  % compute the error at t = T
	  errmax = interpolation_error(@uexact,jmax,npde,t(n),x,u,inf);
	  errL2 = interpolation_error(@uexact,jmax,npde,t(n),x,u);
	  errH1 = interpolation_error(@uexact,jmax,npde,t(n),x,u,2,1,1,@uexactDx);
	  fprintf(' (jmax,errmax,errL2,errH1) = %4d %e %e %e\n', jmax, errmax, errL2, errH1);
	  
   else % for computing error
   
      tspan = [0:0.01:1];
      [t,y]=pdefdm(npde,jmax,tspan,y,yt,@PDE_F,@PDE_G,@BC_L,@BC_R,a,b,reltol,abstol);
      fprintf('\n cpu time used = %e \n', etime(clock,cpu0));
      % compute the error
	  N=size(y,1);
	  Errmax = 0; ErrL2 = 0; ErrH1 = 0;
	  for n=1:N
	     u=reshape(y(n,1:nn),npde,jmax);
	     errmax = interpolation_error(@uexact,jmax,npde,t(n),x,u,inf);
	     errL2 = interpolation_error(@uexact,jmax,npde,t(n),x,u);
	     errH1 = interpolation_error(@uexact,jmax,npde,t(n),x,u,2,1,1,@uexactDx);
	     if n==1
			dt = (t(2)-t(1));
	     elseif n==N
	        dt = (t(N)-t(N-1));
	     else
	        dt = (t(n+1)-t(n-1));
	     end	   
	     Errmax=Errmax+0.5*dt*errmax;
	     ErrL2=ErrL2+0.5*dt*errL2;
	     ErrH1=ErrH1+0.5*dt*errH1;
      end
      fprintf(' (jmax,Errmax,ErrL2,ErrH1) = %4d %e %e %e\n', jmax, Errmax, ErrL2, ErrH1);
	  fprintf(' (jmax,errmax,errL2,errH1) = %4d %e %e %e\n', jmax, errmax, errL2, errH1);
	  
   end

fprintf('cpu time used = %e \n', etime(clock,cpu0));

% -----------------------------------------------------

function f=PDE_F(t,x,u,ux,ut)
  f(1,:) = ut(1,:);
  
function f=PDE_G(t,x,u,ux,ut)
  global epsilon;
  f(1,:) = epsilon*ux(1,:)-u(1,:).*u(1,:)*0.5;
  
function f=BC_L(t,x,u,ux,ut)
%  f(1) = u(1)-uexact(t,x);
  f(1) = ut(1)-uexactDt(t,x);
  
function f=BC_R(t,x,u,ux,ut)
%  f(1) = u(1)-uexact(t,x);
  f(1) = ut(1)-uexactDt(t,x);
    
% -----------------------------------------------------

function u=uexact(t,x)
   global epsilon;
   a1 = (-x+0.5-4.95*t)/(20.0*epsilon);
   a2 = (-x+0.5-0.75*t)/(4.0*epsilon);
   a3 = (-x+0.375)/(2.0*epsilon);
   a = max([a1 a2 a3]);
   r1 = exp((a1-a));
   r2 = exp((a2-a));
   r3 = exp((a3-a));
   u = (0.1*r1+0.5*r2+r3)/(r1+r2+r3);
   
% -----------------------------------------------------

function du=uexactDt(t,x)
   global epsilon;
      a1 = (-x+0.5-4.95*t)/(20.0*epsilon);
      a2 = (-x+0.5-0.75*t)/(4.0*epsilon);
      a3 = (-x+0.375)/(2.0*epsilon);
      a = max([a1, a2, a3]);
      r1 = exp(a1-a);
      r2 = exp(a2-a);
      r3 = exp(a3-a);
      r1t = -r1*4.95/(20*epsilon);
      r2t = -r2*0.75/(4*epsilon);
      r3t = 0.0;
      du = ((r1+r2+r3)*(0.1*r1t+0.5*r2t+r3t)-(r1t+r2t+r3t)*(0.1*r1+0.5*r2+r3)) ...
           /(r1+r2+r3)^2;
      
% -----------------------------------------------------

function du=uexactDx(t,x)
   global epsilon;
      a1 = (-x+0.5-4.95*t)/(20.0*epsilon);
      a2 = (-x+0.5-0.75*t)/(4.0*epsilon);
      a3 = (-x+0.375)/(2.0*epsilon);
      a = max([a1, a2, a3]);
      r1 = exp(a1-a);
      r2 = exp(a2-a);
      r3 = exp(a3-a);
      r1x = -r1/(20*epsilon);
      r2x = -r2/(4*epsilon);
      r3x = -r3/(2*epsilon);
      du = ((r1+r2+r3)*(0.1*r1x+0.5*r2x+r3x)-(r1x+r2x+r3x)*(0.1*r1+0.5*r2+r3)) ...
           /(r1+r2+r3)^2;

