How can I go about solving the following PDE?

Hey, guys

Below is the equation I am trying to solve:

I spent days trying to find appropriate solution adhering to the physical nature of problem. The desired result p(x,t) always give the initial value, and doesn't vary with distance although it should.

Can you give some hints or advice how to go forward? I need fresh advice.

Main file diffusivity_equation.m

*%

 clear all

 clc

 %Parameters shared with the ODE routine

 global n  ncall phi k c cf ct mu h rw q beta xu xl re  B ncase xw IC dx p0 pi

 phi=0.2;

 p0=3000;

 k=125;

 c=1.0e-08;

 B=1.1;

 cf=1.0e-06;

 ct=c+cf;

 mu=1.0;

 h=100;

 xw=3;

 q=10000;

 beta=-(q*B*mu/0.001127*k*h);

 IC=-(q*mu)/(k*h*xw);

 n=41;

 xl=0;

 xu=2000;

 dx=(xu-xl)/(n-1);

 for i=1:n

     x(i)=(i-i)*dx;

 end

 %

 %Independent variable for ODE integration

 t0=0.0;

 tf=12;

 tout=linspace(t0,tf,n);

 nout=n;

 ncall=0;
 %
 p0=3000*ones(n,1);

 %ODE integration

 mf=1;

 reltol=1.0e-04; abstol=1.0e-04;

 options=odeset('RelTol',reltol,'AbsTol',abstol);

 if(mf==1)%FDs

 [t,p]=ode15s(@function_handle2,tout,p0,options);end

 % Display output

 for it=1:nout
    fprintf('\n t = %4.2f\n',t(it));  
 for i=1:n
      fprintf(' x = %4.1f   p(x,t) = %8.5f\n',x(i),p(it,i));
 end

 end
  fprintf('\n ncall = %5d\n',ncall);
 %
 ncall=ncall+1;*

 function_hadle2.m

 function pt=function_handle2(t,p)

 *global  n  ncall phi k c cf ct mu h rw q beta xu xl re  B  px pxx  ndss   ncase 
 xw IC dx

 %

 %Problem parameters

 xl=0.0;

 xu=2000;

 %

 %left boundary condition

 %

 p(1)=0;

 %

 %first-order derivative px

 %

 px=dss004(0,xu,n,p)

 %

 %right boundary condition

 px(n)=0;

 %second-order dervative pxx

 %

 pxx=dss004(0,xu,n,px);

 for i=1:n

    pt(i)=beta*(cf*px(i)^2+pxx(i))

 end 

 %PDE

 pt=pt';

 %

 %Increment calls to function_handle2

 ncall=ncall+1;  *

dss004

 %  File: dss004.m
 %

   function [ux]=dss004(xl,xu,n,u)

   ux(  2)=r4fdx*...
     (  -3.*u(  1) -10.*u(  2) +18.*u(  3)  -6.*u(  4)  +1.*u(  5));

   for i=3:nm2
     ux(  i)=r4fdx*...
     (  +1.*u(i-2)  -8.*u(i-1)  +0.*u(  i)  +8.*u(i+1)  -1.*u(i+2));
   end

   ux(n-1)=r4fdx*...
     (  -1.*u(n-4)  +6.*u(n-3) -18.*u(n-2) +10.*u(n-1)  +3.*u(  n));

   ux(  n)=r4fdx*...
     (   3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u(  n));