% 	pseudo12.m - for Table 10
%This computes mbar(t) 

close('all')

clear all

i=sqrt(-1);

% Set up basic data

% number of points per unit length of (0,a)
hx=1000;

% length of interval
a=10

% total number of points in (0,a)
Nx=hx*a+1;

% set up space variable x

x=zeros(Nx,1);
for r=1:Nx
   x(r,1)=(r-1)/hx;
end;


% specify the initial function f

f=zeros(Nx,1);
for r=1:Nx
   f(r,1)=sin(pi*(r-1)/(Nx-1));
end;
  
  c=20
  %this specifies the amount of accuracy required
  
  b=20
  % this is the diffusion constant
  
  del=0.5-c/(a*b);
  % this is needed to specify the pseudoeigenfunctions
  
  
  %specify the dimension N of the test function space 
  
  M=40  
  N=2*M+1
   
   
% specify the un-normalized pseudo-eigenvectors
      
s=(1:N)-ones(1,N)*(N+1)/2;
   
e=zeros(Nx,N);
for rx=1:Nx              
      e(rx,:)=exp( -0.5*b*x(rx,1)+b*del*x(rx,1)+i*2*pi*s*x(rx,1)/a  )-exp( -0.5*b*x(rx,1)-b*del*x(rx,1)-i*2*pi*s*x(rx,1)/a  );  
end;


%specify the sequence phi
phi=e\f;


%compute  f at time zero
f0=e*phi;



% compute the function ft for all times t from 1 to T

T=10;
mint=zeros(T,1);
maxt=zeros(T,1);
for t=1:T;
   diagt=zeros(N,1);
   for r=1:N
      s=r-(N+1)/2;
      sig=2*pi*s/a;
      diagt(r,1)=exp(t*( -sig^2/b+i*sig-c/a+c^2/(b*a^2)-2*i*sig*c/(a*b)     ));
   end;
   ft=e*(diagt.*phi);
   minf(t,1)=min(real(ft));
   maxf(t,1)=max(real(ft));
end;


% print output

format long
N;
c;

[(1:T)' minf maxf]