%
% solves C_t = (C^mm C_x )_x with given initial condition non-zero only
% on -1 < x < 1
%
% there is similarity solution t^(-1/m+2) F ( x / t^(1/m+2) )
%
%
function nonlindiffusion
%PDEX1  Example 1 for PDEPE
%   This is a simple example with an analytical solution that illustrates 
%   the straightforward formulation, computation, and plotting of the solution
%   of a single PDE.
%
%   In the form expected by PDEPE, the single PDE is
%
%   [pi^2] .*  D_ [u] =  D_ [   Du/Dx   ] +      [0]
%              Dt        Dx
%    ----         ---       -------------    -------------
%     c            u        f(x,t,u,Du/Dx)   s(x,t,u,Du/Dx)
%
%   The equation is to hold on an interval 0 <= x <= 1 for times t >= 0.
%   There is an analytical solution u(x,t) = exp(-t)*sin(pi*x). Boundary
%   conditions are specified so that it will be the solution of the initial-
%   boundary value problem.  Obviously the initial values are sin(pi*x). 
%   Two kinds of boundary conditions are chosen so as to show how they
%   appear in the form expected by the solver.  The left bc is u(0,t) = 0:
%
%     [u]    +     [0] .* [       Du/Dx          ] = [0]
% 
%     ---          ---    ------------------------   ---
%   p(0,t,u)      q(0,t)        f(0,t,u,Du/Dx)        0
%
%   The right bc is taken to be pi*exp(-t) + Du/Dx(1,t) = 0:
%
%   [ pi*exp(-t) ] + [1] .* [   Du/Dx    ] = [0]
%
%   --------------   ---    --------------   ---
%      p(1,t,u)     q(1,t)   f(1,t,u,Du/Dx)   0
%
%   The problem is coded in subfunctions PDEX1PDE, PDEX1IC, and PDEX1BC.
%
%   See also PDEPE, FUNCTION_HANDLE.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 1.6.4.2 $  $Date: 2005/06/21 19:25:04 $

global mm;
m = 0;
mm=input('mm=')
mmbeta=1/(mm+2);
strmmplus2=num2str(mm+2);
strmmpower=['1' '/' strmmplus2];
x = linspace(-20,20,201);
t = linspace(0,40,26);
t1=0:1.6:40;
x1=-10:0.1:10;

nic=input('initial condition? (1 or 2)')

if nic==1
    sol = pdepe(m,@pdex1pde,@pdex1ic1,@pdex1bc,x,t);
else
    sol = pdepe(m,@pdex1pde,@pdex1ic2,@pdex1bc,x,t);
end

        
% Extract the first solution component as u.  This is not necessary
% for a single equation, but makes a point about the form of the output.
u = sol(:,:,1);


% A solution profile can also be illuminating.
figure;
plot(x,u(1,:),'-',x,u(2,:),'-',x,u(10,:),'-',x,u(15,:),'-',x,u(20,:),'-');
axis([-10 10 0 1.1]);
xlabel('x');
ylabel('u(x,t)');
pause;
close;

figure;
plot(x,t1(5)^mmbeta*u(5,:), ...
    '-',x,t1(10)^mmbeta*u(10,:), ...
    '-',x,t1(15)^mmbeta*u(15,:), ...
    '-',x,t1(20)^mmbeta*u(20,:), ... 
    '-',x,t1(25)^mmbeta*u(25,:),'-');
axis([-10 10 0 1]);
xlabel('x');
ylabel('t^{1/m+2}u(x,t)');
pause;
close;

figure;
plot(x./t1(5)^mmbeta, t1(5)^mmbeta*u(5,:),'-', ... 
    x./t1(10)^mmbeta,t1(10)^mmbeta*u(10,:),'-', ... 
    x./t1(15)^mmbeta,t1(15)^mmbeta*u(15,:),'-', ... 
    x./t1(20)^mmbeta,t1(20)^mmbeta*u(20,:),'-', ...
    x./t1(25)^mmbeta,t1(25)^mmbeta*u(25,:),'-');
axis([-10 10 0 1]);
xlabel('x/t^{1/m+2}');
ylabel('t^{1/m+2}u(x/t^{1/m+2},t)');
pause;
close;

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

function [c,f,s] = pdex1pde(x,t,u,DuDx)
global mm;
c = 1;
f = u^mm*DuDx;
s = 0;

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

function u0 = pdex1ic1(x)
if (abs(x)>1)
    u0=0;
else
    u0 = cos(pi*x/2);
end;
    
% --------------------------------------------------------------------------
function u0 = pdex1ic2(x)
if (abs(x-3)<1)
    u0 = cos(pi*(x-3)/2);
else
    if (abs(x+3)<1)
        u0 = cos(pi*(x+3)/2);
    else
        u0=0;
    end
end;
    
% --------------------------------------------------------------------------



function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur;
qr = 0;
  

