function uneven(dt,tf,dx,xu);
% UNEVEN Compares equally- and unequally-spaced finite
% difference methods for the heat equation problem
%    u_t = u_xx, u(0,t)=u(1,t)=0,
%              / 10 sin(33 pi x),   0 < x < 1/3,
%    u(x,0) = <
%              \ sin(3 pi x),   1/3 < x < 1
% Example: 500 steps; nu=0.1 for even; (min nu)=0.4 for uneven
%   >> xu=[0.0:0.01:0.40 0.45:0.05:1.0];  % uneven mesh
%   >> dt=0.00004;  tf=0.02;  dx=0.02;
%   >> uneven(dt,tf,dx,xu);

% ELB 2/25/07

% equally-spaced ("even") method first
x=0:dx:1;  nu=dt/(dx^2);  t=0:dt:tf;
J=length(x)-1;  % J is number of space steps for even
UE=zeros(size(x));  % initial condition in pieces:
UE(x<=1/3)=10*sin(33*pi*x(x<=1/3));   UE(x>1/3)=sin(3*pi*x(x>1/3));
erre(1)=0;
for n=1:length(t)-1
    UE(2:J)=UE(2:J) + nu*( UE(3:J+1) - 2*UE(2:J) + UE(1:J-1) );
    erre(n+1)=max(abs(exactsoln(x,n*dt)-UE));
end

% Unequally-spaced method next
JU=length(xu)-1;  % xu is length JU+1
dxU=xu(2:JU+1)-xu(1:JU);  % dxU is length JU
UU=zeros(size(xu));
UU(xu<=1/3)=10*sin(33*pi*xu(xu<=1/3));  UU(xu>1/3)=sin(3*pi*xu(xu>1/3));
erru(1)=0;
for n=1:length(t)-1
    nuU=2*dt./(dxU(1:JU-1)+dxU(2:JU));
    dRight=(UU(3:JU+1) - UU(2:JU))./dxU(2:JU);
    dLeft=(UU(2:JU) - UU(1:JU-1))./dxU(1:JU-1);
    UU(2:JU)=UU(2:JU) + nuU .* (dRight - dLeft);
    erru(n+1)=max(abs(exactsoln(xu,n*dt)-UU));
end

% plot solutions; plot errors
JJ=2000;  xx=linspace(0,1,JJ);
figure(1), clf, plot(x,UE,'+',xu,UU,'o',xx,exactsoln(xx,tf))
xlabel x, ylabel u, legend('equally spaced','unequally-spaced','exact')
title(['exact and approx solns at t_f = ' num2str(tf)])
figure(2), clf, semilogy(t,erre,'+',t,erru,'o')
xlabel t, ylabel('error'), legend('equally spaced','unequally-spaced')

    function u=exactsoln(x,t);
        % subfunction to compute the exact soln
        N=500; m=1:N; a=zeros(1,N);
        a(3)=2/3; a(33)=10/3;
        mm=m((m~=3)&(m~=33));  % indices for which general formula applies
        temp=10*(sin((33-mm)*pi/3)./((33-mm)*pi)-sin((33+mm)*pi/3)./((33+mm)*pi));
        a(mm)=temp+sin((3+mm)*pi/3)./((3+mm)*pi) - sin((3-mm)*pi/3)./((3-mm)*pi);
        pim=pi*(1:N);  ae=a.*exp(-t*(pim.^2));
        u=sum(repmat(ae',1,length(x)).*sin(pim'*x));
    end
end