function bkfplt = bkfplt(x,j,Theta,v,color)

%
%       bkfplt = bkfplt(x,j,Theta[,v[,color]])
%
%  plot the (Birkhoff) spline which represents the j-th derivative at x of the
%  error in Hf the Hermite interpolant to f at Theta a multiset of n points, as:
%
%       D^j f(x) - D^j Hf(x) = integral spline(t)*D^n f(t) dt
%  
%  the domain of the plot is the convex hull of {x,Theta} (which is the support %  of the spline if it is non-zero) unless it is specified by the optional 
%  4-th argument v=[xmin xmax]. By specifying v it is possible to zoom up on 
%  places where the spline behaves in an interesting way
%
%  The points in Theta are plotted as 'o' (with those repeated placed on top of
%  each other), x is a ',' and the Chebshev intervals are plotted as lines
%
%  The optional 5-th argument if present causes the spline to be plotted in 
%  colour:  white plot, yellow Theta, red x, green Chebyshev intervals
%
%  example:  spline = bkfplt(0.67,2,[0 0.2 0.4 0.6 0.8 1])
%
%  See also:  bkfspmak (which constructs the spline which is plotted)
%
%  S. Waldron October 1994
%

n=length(Theta);   % number of points in Theta

% limits of the co-ordintate axes
if nargin>=4
   xmin=v(1);   xmax=v(2);
else
   xmin=min([Theta,x]);   xmax=max([Theta,x]);
end

% decide whether to plot in colour
if nargin>4
   plotcolor='w'; thetacolor='y'; xcolor='r'; intervalcolor='g';
   % white plot, yellow Theta, red x, green Chebyshev intervals
else
   plotcolor='w'; thetacolor='w'; xcolor='w'; intervalcolor='w';
   % everything is plotted as white
end
   
% make the spline
M=bkfspmak(x,j,Theta);

% clear graphics screen
clf reset

% plot (in white) 101 points from xmin to xmax
h=(xmax-xmin)/100;
xx=[xmin:h:xmax];  yy=fnval(M,xx);
plot(xx,fnval(M,xx),'color',plotcolor)
hold on
ymin=min(yy);  ymax=max(yy);

pause

% display the points in Theta
if ymax-ymin>0
   ylength=ymax-ymin;
else 
   % the spline is zero, and we plot the lower and uppermost endpoints of the
   % the current plot in black and make the current y range the ylength of the
   % spline
   v=axis;  
   ylength=v(4)-v(3);
   plot(v(1),v(3),'color','k','linestyle','.')
   plot(v(2),v(4),'color','k','linestyle','.')
end
% plot the points in Theta as (yellow) 'o'
% with points that occur with multiplicity r placed on top of each other
% (so that a pile of 35 would have the same length as the ylength of the spline)
hh=ylength/35;  % how high to plot the next copy of a point above the last
for i=1:n
   theta=Theta(i);
   if (theta>=xmin) & (theta<=xmax)  % if theta from the domain we are plotting
      % find r the multiplicity of theta in Theta
      r=n-length(find(Theta-theta*ones(size(Theta)))); 
      % plot r copies of theta
      for k=1:r
         plot(theta,(k-1)*hh,'linestyle','o','color',thetacolor,'markersize',5)
      end
   end
end

% plot x
plot(x,0,'linestyle','.','color',xcolor,'markersize',19)

print -deps2 aspline
%print -depsc2 aspline
%print -dps2 aspline
pause

% display J:=J(Theta) the union of n-j disjoint open intervals with total length
% j/(n-1)*diam(Theta) contained within conv(Theta) mentioned in the paper of
% Shadrin. He predicts the Birkhoff spline has exactly one zero in the interior
% of conv(Theta) when x is in J, and otherwise is of constant sign
[a,b]=bkfchsg(x,j,Theta);   % find the endpoints of the n-j intervals
for i=1:n-j
   plot([a(i) b(i)],[0 0],'color',intervalcolor)  
   %plot([a(i) b(i)],[i/10 i/10],'color',intervalcolor) 
   % plot the intervals at heights 0.1,0.2,...
end