function rd1d

    name = [];
    while (isempty(name))
        clc
        disp('Reaction-Diffusion Lab (1d)')
        disp('---------------------------')
        disp('(Version Date: 30 Nov 2001)')
        disp(' ')
        disp('First enter your name.  This will allow you to identify any plots that you print.')
        name = input('Name: ','s');
    end
    while (1)
        c = [];
        while (isempty(c))
            clc
            disp('Reaction-Diffusion Lab (1d)')
            disp('---------------------------')
            disp(' ')
            str = sprintf('Welcome, %s',name);
            disp(str)
            disp(' ')
            disp('Options are:')
            disp('1.  Diffusion only')
            disp('2.  Reaction only ("well mixed")')
            disp('3.  Both reaction and diffusion')
%            disp('4.  Information about the equations and this program')
            disp(' ')
            c = input('Enter choice (1, 2, 3, or q to quit): ','s');
        end
        if (c == 'q')
            clc
            return
        elseif (c == '1')
            rd1d_diffusion(name)
        elseif (c == '2')
            rd1d_reaction(name)
        elseif (c == '3')
            rd1d_both(name)
            %        elseif (c == '4')
            %rd1d_info
        end
    end


function rd1d_info
    clc
    disp('INFORMATION about the equations and this program')
    disp('------------------------------------------------')
    disp(' ')
    disp('Not implemented yet!')
    input('\nHit Enter to continue.','s')
    return

function rd1d_diffusion(name)

    while (1)
        c = [];
        while (isempty(c))
            clc
            disp('DIFFUSION ONLY')
            disp('--------------')
            disp('The goal here is simply to observe the diffusion process,')
            disp('and see how the behavior of a diffusing chemical depends')
            disp('on the diffusion coefficient.')
            disp(' ')
            disp('There are two demonstrations:')
            disp('1.  Diffusion of random initial conditions.')
            disp('2.  Comparison of different diffusion coefficients.')
            c = input('Enter choice (1, 2, or q to quit): ','s');
        end
        if (c == 'q')
            clc
            return
        elseif (c == '1')
            rd1d_diffusion_difftimes(name)
        elseif (c == '2')
            rd1d_diffusion_diffcoeffs(name)
        end
    end

    
function rd1d_diffusion_diffcoeffs(name)

    N = 100;
    L = 5;
    x = linspace(0,L,N+1)';
    h = L/N;
    T = 1;


    a = 0;
    b = 0;
    g = 0;
    d2 = 0;

    u0 = 0.2*ones(size(x));
    u0(49) = 2.4;
    u0(50) = 2.4;
    u0(51) = 2.4;
    v0 = zeros(size(x));


    z0 = [u0; v0];
    
    disp('Working...')

    opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
    ts = linspace(0,T,101);
    [t1,z1] = ode15s(@rdsys0,ts,z0,opts,[0 0 0 .01 0 h]);
    [t2,z2] = ode15s(@rdsys0,ts,z0,opts,[0 0 0 .1 0 h]);
    [t3,z3] = ode15s(@rdsys0,ts,z0,opts,[0 0 0 1 0 h]);
    
    disp(' ')
    input('Hit Enter to see how the concentrations change over time.','s');
    
    figure(1)
    clf
    set(gcf,'DoubleBuffer','on');
    lw = 2.5;
    axmax = 2.5; % for now...

    M = length(t1);
    
    for k = 1:M,
        hold off
        plot(x,u0,'b:')
        hold on
        plot(x,z1(k,1:N+1),'b','Linewidth',[lw]);
        plot(x,z2(k,1:N+1),'b--','Linewidth',[lw]);
        plot(x,z3(k,1:N+1),'b-.','Linewidth',[lw]);
        axis([ 0 L 0 axmax])
        str = sprintf('t=%9.5f',t1(k));
        title(str)
        xlabel('x')
        ylabel('concentration')
        drawnow;
    end

    str = sprintf('Name: %s',name);
    text(3,0.1,str)
    title('Diffusion Demonstration #2')
    text(0.1,2.35,'Same starting concentrations (thin dashed line).');
    text(0.1,2.20,'Same ending time (t=1).');
    text(0.1,2.05,'Three different diffusion coefficients: 0.01, 0.1, and 1.0.');
    text(0.1,1.90,'Which is which?');

    input('\nHit Enter to continue.','s');
    return

    
function rd1d_diffusion_difftimes(name)

    N = 100;
    L = 5;
    x = linspace(0,L,N+1)';
    h = L/N;
    
    a = 0;
    b = 0;
    g = 0;
    d1 = 1;
    d2 = 0;

    u0 = 0.75+x/10+randpert(0.2,length(x));
    v0 = zeros(size(x));

    figure(1)
    clf
    set(gcf,'DoubleBuffer','on');
    lw = 2.5;
    axmax = 2; % for now...

    z0 = [u0; v0];
    
    disp('Working...')

    opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
    T1 = 0.2;
    T2 = 0.6;
    T3 = 1.8;
    T4 = 5.4;
    [t1,z1] = ode15s(@rdsys0,[0 T1],z0,opts,[0 0 0 .01 0 h]);
    [t2,z2] = ode15s(@rdsys0,[0 T2],z0,opts,[0 0 0 .1 0 h]);
    [t3,z3] = ode15s(@rdsys0,[0 T3],z0,opts,[0 0 0 1 0 h]);
    [t4,z4] = ode15s(@rdsys0,[0 T4],z0,opts,[0 0 0 1 0 h]);
    
    hold off
    M1 = length(t1);
    styles = {'b:','b','b-.','b--'};
    rstyles = styles(randperm(4));
    plot(x,z1(M1,1:N+1),char(rstyles(1)),'Linewidth',[lw]);
    hold on
    M2 = length(t2);
    plot(x,z2(M2,1:N+1),char(rstyles(2)),'Linewidth',[lw]);
    M3 = length(t3);
    plot(x,z3(M3,1:N+1),char(rstyles(3)),'Linewidth',[lw]);
    M4 = length(t4);
    plot(x,z4(M4,1:N+1),char(rstyles(4)),'Linewidth',[lw]);
    axis([ 0 L 0 axmax])
    xlabel('x')
    ylabel('concentration')
    str = sprintf('Name: %s',name);
    text(3,0.1,str)
    title('Diffusion Demonstration #1')
    text(0.1,1.90,'Same starting concentration (not shown).');
    text(0.1,1.75,'Same diffusion coefficient.');
    tstr = sprintf('Four different ending times: t=%3.1f, t=%3.1f, t=%3.1f and t=%3.1f.',T1,T2,T3,T4);
    text(0.1,1.60,tstr);
    text(0.1,1.45,'Which is which?');

    input('\nHit Enter to continue.','s');
    return

    
function rd1d_reaction(name)

    clc
    disp('REACTION ONLY')
    disp('-------------')
    disp(' ')
    disp('If the reactants are "well mixed", then their concentrations')
    disp('are uniform.  That is, they are independent of x.  In this case,')
    disp('their is no diffusion, and the behavior of the system is')
    disp('governed by the equations')
    disp(' ')
    disp('  du/dt = g*(a - u + u^2*v)')
    disp('  dv/dt = g*(b - u^2*v)')
    disp(' ')
    disp('We will use a = 0.025, b = 1.1, and g = 4.')
    disp(' ')
    disp('Let''s run a simulation of this system...')
    disp(' ')
    rd1d_reaction_sim(name)
    input('\nHit Enter to continue.','s');
    return

    
function rd1d_reaction_sim(name)

    a = 0.025;
    b = 1.1;
    g = 4;
    p = [a b g];
    
    u0 = inputnumber('Enter the initial value of u (between 0 and 3, or q to quit): ',[0 3],[1 1]);
    if (isnan(u0))
        return
    end
    v0 = inputnumber('Enter the initial value of v (between 0 and 3, or q to quit): ',[0 3],[1 1]);
    if (isnan(v0))
        return
    end

%    T = inputnumber('Enter the time duration for the simulation (between 0 and 20, or q to quit): ', ...
%                            [0 20],[0 1]);
    T = 10;
    disp('Working...')
    [t,z] = ode15s(@reactionodes,[0 T],[u0 v0]',[],p);
    figure(1)
    clf
    set(gcf,'DoubleBuffer','on')
    subplot(2,1,1)
    plot(t,z(:,1),'b',t,z(:,2),'k--','Linewidth',[2.5])
    xlabel('t')
    ylabel('concentration')
    legend('u','v')
    title('Reaction Only: u and v vs. t')
    hold on
    plot(t,(b/(a+b)^2)*ones(size(t)),'k--');
    plot(t,(b+a)*ones(size(t)),'b');
    hold off
    ax1 = axis;
    
    subplot(2,1,1)
    h = line([0 0],[ax1(3) ax1(4)],'Color',[.5 .5 .5]);
    M = length(t);
    for k = 1:M,
        subplot(2,1,2)
        plot([0 5],[z(k,1) z(k,1)],'b',[0 5],[z(k,2) z(k,2)],'k--','Linewidth',[2.5])
        axis([0 5 0 ax1(4)])
        xlabel('x')
        ylabel('concentration')
        str = sprintf('t=%9.5f',t(k));
        title(str)
        set(h,'XData',[t(k) t(k)])
        pause(0.05)
        drawnow
    end
    delete(h)
    str = sprintf('Reaction Only: Concentrations u and v vs. x at t=%4.1f',T);
    title(str)
    str = sprintf('Name: %s',name);
    text(3,0.15,str)

    return
    
    
function rd1d_both(name)

    while (1)
        c = [];
        while (isempty(c))
            clc
            disp('REACTION-DIFFUSION')
            disp('------------------')
            disp('Here we include both the reaction and diffusion processes.')
            disp('The goal will be to investigate two questions:')
            disp('1.   How does changing the diffusion coefficient of one of')
            disp('     the reactants affect the long-term behavior of the system?')
            disp('2.   How does changing the length of the interval change')
            disp('     the long-term behavior of the system?')
            c = input('\nChoose one of the questions (enter 1, 2, or q to quit): ','s');
        end
        if (c == 'q')
            return
        elseif (c == '1')
            rd1d_both_diffcoeff(name)
        elseif (c == '2')
            rd1d_both_length(name)
        end
    end

    
function rd1d_both_diffcoeff(name)

    disp(' ')
    disp('The diffusion coefficient of u is 1.');
    disp('Here you can choose the diffusion coefficient for v.');
    d2 = inputnumber('Enter the diffusion coefficient for v (a number between 0 and 100, or q to quit): ',[0 100],[0 1]);
    if (isnan(d2))
        return
    end
    
    N = 100;
    L = 5;
    x = linspace(0,L,N+1)';
    h = L/N;
    Tmax = 500;


    a = 0.025;
    b = 1.1;
    g = 4;
    d1 = 1;

    u0 = (a+b)*ones(size(x)) + randpert(0.1,length(x));
    v0 = b/(a+b)^2 * ones(size(x)) + randpert(0.1,length(x));
    figure(1)
    clf
    set(gcf,'DoubleBuffer','on');
    lw = 2.5;
    axmax = 3; % for now...

    hold off
    hv = plot(x,v0,'k--','Linewidth',[lw]);
    hold on
    plot(x,(b/(a+b)^2)*ones(size(x)),'k--');
    hu = plot(x,u0,'b','Linewidth',[lw]);
    plot(x,(b+a)*ones(size(x)),'b');
    axis([ 0 L 0 axmax])
    xlabel('x')
    ylabel('concentration')

    legend([hu hv],'u','v')
    title('Initial Concentrations for the Reaction-Diffusion System')
    str = sprintf('Diffusion Coefficient of v is %5.2f',d2);
    text(0.1,0.95*axmax,str)
    str = sprintf('Name: %s',name);
    text(3,0.1,str)
    
    disp(' ')
    disp('Figure No. 1 shows the initial random concentrations of the reactants.');
    disp('The thin lines show the equilibrium concentrations for the "well-mixed" system.');
    
    disp('Move the cursor into the graphics window, and use');
    disp('the mouse to modify the initial concentrations.');
    disp('Use the left button for u, and the right button for v.');
    disp('Hit Enter to begin computing the solution.');

    while (1)
        [x0,y0,btn] = ginput(1);
        if (isempty(x0))
            break
        end
        if (x0 < 0 | x0 > L | y0 < 0 | y0 > axmax)
            continue;
        end
        if (btn == 1)
            u1 = spline(x,[0 u0' 0],x0);
            s = (y0-u1)/3;
            r = 1+s*sech(x-x0).^2;
            u0 = u0.*r;
        else
            v1 = spline(x,[0 v0' 0],x0);
            s = (y0-v1)/3;
            r = 1+s*sech(x-x0).^2;
            v0 = v0.*r;
        end
        hold off
        hv = plot(x,v0,'k--','Linewidth',[lw]);
        hold on
        plot(x,(b/(a+b)^2)*ones(size(x)),'k--');
        hu = plot(x,u0,'b','Linewidth',[lw]);
        plot(x,(b+a)*ones(size(x)),'b');
        axis([ 0 L 0 axmax])
        xlabel('x')
        ylabel('concentration')

        legend([hu hv],'u','v')
        title('Initial Concentrations for the Reaction-Diffusion System')
        str = sprintf('Diffusion Coefficient of v is %5.2f',d2);
        text(0.1,0.95*axmax,str)
        str = sprintf('Name: %s',name);
        text(3,0.1,str)

    end
 
    z0 = [u0; v0];
    
    p = [a b g d1 d2 h]';

    disp('Working...')

    [t,z] = solveit(z0,p,5,Tmax,1e-5);
    
    mm = max(max(z));
    if (mm < 3)
        axmax = 3;
    else
        axmax = 1.1*mm;
    end

    M = length(t);
    
    disp(' ')
    input('Hit Enter to see how the concentrations change over time.','s');
    
    for k = 2:M,
        hold off
        hv = plot(x,z(k,N+2:2*(N+1)),'k--','Linewidth',[lw]);
        hold on
        plot(x,(b/(a+b)^2)*ones(size(x)),'k--');
        hu = plot(x,z(k,1:N+1),'b','Linewidth',[lw]);
        plot(x,(b+a)*ones(size(x)),'b');
        axis([ 0 L 0 axmax])
        str = sprintf('t=%9.5f',t(k));
        title(str)
        xlabel('x')
        ylabel('concentration')
        drawnow;
    end

    legend([hu hv],'u','v')
    title('Equilibrium Solution of the Reaction-Diffusion System')
    str = sprintf('Diffusion Coefficient of v is %5.2f',d2);
    text(0.1,0.95*axmax,str)
    str = sprintf('Name: %s',name);
    text(3,0.1,str)
    disp('Done.')
    return


function rd1d_both_length(name)

    disp(' ')
    disp('The diffusion coefficient of u is 1, and the diffusion coefficient of v is 10.');
    disp('Here you can choose the length of the domain, L.');
    L = inputnumber('Enter the domain length L ( 0 < L <= 40, or q to quit): ',[0 40],[0 1]);
    if (isnan(L))
        return
    end
    
    N = max(ceil(L)*5,24);
    x = linspace(0,L,N+1)';
    h = L/N;
    Tmax = 500;

    a = 0.025;
    b = 1.1;
    g = 4;
    d1 = 1;
    d2 = 10;

    u0 = (a+b)*ones(size(x)) + randpert(0.1,length(x));
    v0 = b/(a+b)^2 * ones(size(x)) + randpert(0.1,length(x));

    figure(1)
    clf
    set(gcf,'DoubleBuffer','on');
    lw = 2.5;
    axmax = 3; % for now...

    hold off
    hv = plot(x,v0,'k--','Linewidth',[lw]);
    hold on
    plot(x,(b/(a+b)^2)*ones(size(x)),'k--');
    hu = plot(x,u0,'b','Linewidth',[lw]);
    plot(x,(b+a)*ones(size(x)),'b');
    axis([ 0 L 0 axmax])
    xlabel('x')
    ylabel('concentration')

    legend([hu hv],'u','v')
    title('Initial Concentrations for the Reaction-Diffusion System')
    str = sprintf('Length of the domain is %6.2f',L);
    text(0.02*L,0.95*axmax,str)

    str = sprintf('Name: %s',name);
    text(0.6*L,0.1,str)
    
    disp(' ')
    disp('Figure No. 1 shows the initial random concentrations of the reactants.');
    disp('The thin lines show the equilibrium concentrations for the "well-mixed" system.');
    
    disp('Move the cursor into the graphics window, and use');
    disp('the mouse to modify the initial concentrations.');
    disp('Use the left button for u, and the right button for v.');
    disp('Hit Enter to begin computing the solution.');

    while (1)
        [x0,y0,btn] = ginput(1);
        if (isempty(x0))
            break
        end
        if (x0 < 0 | x0 > L | y0 < 0 | y0 > axmax)
            continue;
        end
        if (btn == 1)
            u1 = spline(x,[0 u0' 0],x0);
            s = (y0-u1)/3;
            r = 1+s*sech((x-x0)*5/L).^2;
            u0 = u0.*r;
        else
            v1 = spline(x,[0 v0' 0],x0);
            s = (y0-v1)/3;
            r = 1+s*sech((x-x0)*5/L).^2;
            v0 = v0.*r;
        end
        hold off
        hv = plot(x,v0,'k--','Linewidth',[lw]);
        hold on
        plot(x,(b/(a+b)^2)*ones(size(x)),'k--');
        hu = plot(x,u0,'b','Linewidth',[lw]);
        plot(x,(b+a)*ones(size(x)),'b');
        axis([ 0 L 0 axmax])
        xlabel('x')
        ylabel('concentration')

        legend([hu hv],'u','v')
        title('Initial Concentrations for the Reaction-Diffusion System')
        str = sprintf('Length of the domain is %6.2f',L);
        text(0.02*L,0.95*axmax,str)
        str = sprintf('Name: %s',name);
        text(0.6*L,0.1,str)

    end
 
    z0 = [u0; v0];
    
    p = [a b g d1 d2 h]';

    disp('Working...')

    [t,z] = solveit(z0,p,5,Tmax,1e-5);
    
    mm = max(max(z));
    if (mm < 3)
        axmax = 3;
    else
        axmax = 1.1*mm;
    end

    M = length(t);
    
    disp(' ')
    input('Hit Enter to see how the concentrations change over time.','s');
    
    for k = 2:M,
        hold off
        hv = plot(x,z(k,N+2:2*(N+1)),'k--','Linewidth',[lw]);
        hold on
        plot(x,(b/(a+b)^2)*ones(size(x)),'k--');
        hu = plot(x,z(k,1:N+1),'b','Linewidth',[lw]);
        plot(x,(b+a)*ones(size(x)),'b');
        axis([ 0 L 0 axmax])
        str = sprintf('t=%9.5f',t(k));
        title(str)
        xlabel('x')
        ylabel('concentration')
        drawnow;
    end

    legend([hu hv],'u','v')
    title('Equilibrium Solution of the Reaction-Diffusion System')
    str = sprintf('Length of the domain is %6.2f',L);
    text(0.02*L,0.95*axmax,str)
    str = sprintf('Name: %s',name);
    text(0.6*L,0.1,str)
    disp('Done.')
    return

    
function [t,z] = solveit(z0,p,Tstep,Tmax,eqtol)

    format compact
    t = [];
    z = [];
    opts = odeset('RelTol',1e-5,'AbsTol',1e-6);
    [t,z] = ode15s(@rdsys0,[0 Tstep],z0,opts,p);
    disp(Tstep);
    m = length(t);
    done = isequil(z(m,:)',p,eqtol);
    while (~done)
        [t1,z1] = ode15s(@rdsys0,linspace(t(m),t(m)+Tstep,21),z(m,:),opts,p);
        disp(t(m)+Tstep);
        t = [t; t1];
        z = [z; z1];
        m = length(t);
        done = isequil(z(m,:)',p,eqtol) | t(length(t)) >= Tmax;
    end
    if (t(length(t)) >= Tmax)
        wstr = sprintf('WARNING: The solution has not converged to an equilibrium after t=%7.2f',Tmax);
        disp(wstr)
    end

    
function q = isequil(z0,p,tol)

    e = max(abs(rdsys0(0,z0,p)));
    q = e < tol;

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

function w = inputnumber(promptstr, bounds, equalallowed)

    wmin = bounds(1);
    wmax = bounds(2);
    wstr = input(promptstr,'s');
    w = str2double(wstr);
    while (isempty(wstr) | (isnan(w) & wstr ~= 'q' ) | w < wmin | w > wmax | ...
                         (w==wmin & ~equalallowed(1)) | (w==wmax & ~equalallowed(2)))
        beep;
        disp('Sorry, that is not a valid number.')
        wstr = input(promptstr,'s');
        w = str2double(wstr);
    end
    return

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

function p = randpert(a,n)
   p1 = a*cumsum(rand(n,1)-0.5);
   p = p1 - mean(p1);
   return
   

function dz = rdsys0(t,z,p)
   a = p(1);
   b = p(2);
   g = p(3);
   d1 = p(4);
   d2 = p(5);
   h = p(6);
   N = length(z)/2;
   u = z(1:N);
   v = z(N+1:2*N);
   du = g*(a-u+u.^2.*v) + d1*diff2(u,h,0);
   dv = g*(b-u.^2.*v)   + d2*diff2(v,h,0);
   dz = [du; dv];
   return
   

function d2h = diff2(x,h,bcs)
   N = length(x);
   if (bcs == 0)
      xu = [x(1); x(1:N-1)];
      xd = [x(2:N); x(N)];
   else
      xu = [x(N); x(1:N-1)];
      xd = [x(2:N); x(1)];
   end
   d2h = (xd - 2*x + xu)/h^2;
   return
   

function dz = reactionodes(t,z,p)
   a = p(1);
   b = p(2);
   g = p(3);
   u = z(1);
   v = z(2);
   du = g*(a - u + u^2*v);
   dv = g*(b - u^2*v);
   dz = [du; dv];
   return