%
% Clear previous files
  clear all
  clc
%
% Parameters shared with the ODE routine
  global xl xu nx ncall ndss ncase
%
% Grid in x, initial condition
  xl=0;xu=1;nx=26;
  for ix=1:nx
    xg(ix)=(xu-xl)*(ix-1)/(nx-1);
    u0(ix)=0; 
  end   
%
% Independent variable for ODEs
  t0=0.0;tf=1.5;nout=11;
  tout=linspace(t0,tf,nout); 
  ncall=0;
%
% ODE integration
  ncase=1;ndss=4;
  [t,u]=ode15s(@pde_1,tout,u0);
%
% Display selected output
  for it=1:nout
    if(it==1);u(it,1)=0;end
    if(it> 1);u(it,1)=1;end
    fprintf('\n      t       x      u(x,t)\n');
    for ix=1:5:nx
      fprintf('%7.2f%8.2f%12.5f\n',t(it),xg(ix),u(it,ix));
    end
  end  
  fprintf('\n ncall = %4d\n',ncall);
%
% Plot numerical solution
  figure(1);
  plot(xg,u);
  xlabel('x');ylabel('u(x,t)');
  title('u(x,t) vs x, t=0,0.15,...,1.5');
