How to stabilize discretized ODEs by method of lines?

Hello, guys

I am about to solve linear diffusivity equation representing pressure distribution in reservoir for transient conditions. Originally, the pressure should gradually i.e. exponentially rise as it goes further from the center of the well. However, matlab shows pressure declining tendency as it goes from the well. How can I stabilize the solution so that I can obtain proper solution? The correct solution should be pressure gradually increasing from x=0 to x=1000.

The codes are attached below:

 *% lineardiffusion1.m*

 global ncall D x0 IP  

 D=3000;
 x0=1000;
 n=51;
 dx=x0/(n-1);
 IP=3500;
 u0=IP*ones(n,1);
 for i=1:n
    x(i)=(i-1)*dx+50;
end
tf=50
tout=[0.0:2:tf]';
nout=length(tout);
ncall=0;
reltol=1.0e-06; abstol=1.0e-06;
options=odeset('RelTol',reltol,'AbsTol',abstol);
[t,u]=ode15s(@lineardiffusion,tout,u0,options)

 %
 % Display a heading and selected output
  fprintf('\n nx = %2d   x0 = %4.2f   std = %4.2f\n',nx,x0,std);
  disp('at one')
 % pause
  for it=1:nout
      it
    fprintf('\n t = %4.2f\n',t(it));  
    disp('at two')
  %pause
    for i=1:2:nx
      fprintf(' x = %4.1f   u(x,t) = %8.5f\n',x(i),u(it,i));
    end
  end
  fprintf('\n ncall = %5d\n',ncall);
%
% Parametric plot
  figure(1);
  plot(x,u); axis([0 x0 0 10000]);
  title('u(x,t), t = 0, 0.25, ..., 2.5'); xlabel('x'); ylabel('u(x,t)')

surf(x, t ,u);
%contour(r, tt, u2);
axis([0 x0 0 tf 0 10000]);
xlabel('Linear distance');
ylabel('time');
zlabel('pressure');
title([{'Surface plot at time t=',tf}])%, ...
                     % { 'at time t=7.5'}]);
rotate3d on
disp('check plot and press return to progress')
pause
hold 'off'

% Create figure
figure2 = figure('PaperSize',[20.98 29.68]);

 % Create axes
axes2 = axes('Parent',figure2);
box('on');
hold('all');

si=size(u,1);
plot(axes2,x,u(1,:),'k.')
time(1)=0;
hold on
plot(axes2,x,u(10,:),'m.')
time(2)=t(10);
plot(x,u(15,:),'b.')
time(3)=t(15);
plot(x,u(20,:),'g.')
time(4)=t(20);
plot(x,u(30,:),'c.')
time(4)=t(30);
plot(x,u(si,:),'r.')
time(5)=t(si);

 % Create legend

 legend2=legend(axes2,'t=0','t=0.225','t=0.625','t=1.8725','t=2.5');    % this works

 set(legend2,'Position',[0.5354 0.2781 0.1946 0.2333]);

 print -dpdf 'doc.pdf'

 *%lineardiffusion*

function ut=lineardiffusion(t,u)
global ncall D x0 IP xl xu nx
nx=51;
for i=1:nx
  ux=dss008(0.0,x0,nx,u);
  u(nx)=IP;   
  uxx=dss008(0.0,x0,nx,ux);
  ux(1)=-100;
  ut(i)=D*(uxx(i));
end
ut=ut';
ut(nx)=0;
ncall=ncall+1;

 %  File: dss008.m
 %
    function [ux]=dss008(xl,xu,n,u)

 %  Compute the spatial increment
    dx=(xu-xl)/(n-1);
    r8fdx=1./(40320.*dx);
    nm4=n-4;
 %
 %  Grid point 1
    ux(  1)=r8fdx*...
     (-109584.     *u(  1)...
      +322560.     *u(  2)...
      -564480.     *u(  3)...
      +752640.     *u(  4)...
      -705600.     *u(  5)...
      +451584.     *u(  6)...
      -188160.     *u(  7)...
      +46080.      *u(  8)...
      -5040.       *u(  9));
%
%  Grid point 2
   ux(  2)=r8fdx*...
     (-5040.       *u(  1)...
      -64224.      *u(  2)...
      +141120.     *u(  3)...
      -141120.     *u(  4)...
      +117600.     *u(  5)...
      -70560.      *u(  6)...
      +28224.      *u(  7)...
      -6720.       *u(  8)...
      +720.        *u(  9));
%
%  Grid point 3
   ux(  3)=r8fdx*...
     (+720.        *u(  1)...
      -11520.      *u(  2)...
      -38304.      *u(  3)...
      +80640.      *u(  4)...
      -50400.      *u(  5)...
      +26880.      *u(  6)...
      -10080.      *u(  7)...
      +2304.       *u(  8)...
      -240.        *u(  9));
%
%  Grid point 4
   ux(  4)=r8fdx*...
     (-240.        *u(  1)...
      +2880.       *u(  2)...
      -20160.      *u(  3)...
      -18144.      *u(  4)...
      +50400.      *u(  5)...
      -20160.      *u(  6)...
      +6720.       *u(  7)...
      -1440.       *u(  8)...
      +144.        *u(  9));
%
%  Grid point i
   for i=5:nm4
     ux(  i)=r8fdx*...
     (+144.        *u(i-4)...
      -1536.       *u(i-3)...
      +8064.       *u(i-2)...
      -32256.      *u(i-1)...
      +0.          *u(i  )...
      +32256.      *u(i+1)...
      -8064.       *u(i+2)...
      +1536.       *u(i+3)...
      -144.        *u(i+4));
   end
%
%  Grid point n-3
   ux(n-3)=r8fdx*...
     (-144.        *u(n-8)...
      +1440.       *u(n-7)...
      -6720.       *u(n-6)...
      +20160.      *u(n-5)...
      -50400.      *u(n-4)...
      +18144.      *u(n-3)...
      +20160.      *u(n-2)...
      -2880.       *u(n-1)...
      +240.        *u(n  ));
%
%  Grid point n-2
   ux(n-2)=r8fdx*...
     (+240.        *u(n-8)...
      -2304.       *u(n-7)...
      +10080.      *u(n-6)...
      -26880.      *u(n-5)...
      +50400.      *u(n-4)...
      -80640.      *u(n-3)...
      +38304.      *u(n-2)...
      +11520.      *u(n-1)...
      -720.        *u(n  ));
%
%  Grid point n-1
   ux(n-1)=r8fdx*...
     (-720.        *u(n-8)...
      +6720.       *u(n-7)...
      -28224.      *u(n-6)...
      +70560.      *u(n-5)...
      -117600.     *u(n-4)...
      +141120.     *u(n-3)...
      -141120.     *u(n-2)...
      +64224.      *u(n-1)...
      +5040.       *u(n  ));
%
%  Grid point n
   ux(n  )=r8fdx*...
     (+5040.       *u(n-8)...
      -46080.      *u(n-7)...
      +188160.     *u(n-6)...
      -451584.     *u(n-5)...
      +705600.     *u(n-4)...
      -752640.     *u(n-3)...
      +564480.     *u(n-2)...
      -322560.     *u(n-1)...
      +109584.     *u(n  ));