% EXAMPLE3_55 Computations for Example 3.55.
%
f1='%10.0f %12.8f %12.8f\n';
N=40;
load -ascii abhrhermite;
ab0=abhrhermite(1:N+2,:); 
ab1=chri1(N+1,ab0,0);
ab=chri1(N,ab1,0); ab(1,2)=sqrt(pi)/2;
n=5;
xw=gauss(n,ab);
xw(:,2)=xw(:,2)./xw(:,1);
for k=1:n
  fprintf(f1,k,xw(k,1),xw(k,2))
end
fprintf('\n')
nx=500;
x=linspace(0,max(xw(:,1)),nx);
y=zeros(1,nx);
y=stepfunc(n,x,xw);
hold on
ymaxwell=exp(-x.^2);
plot(x',ymaxwell')
plot(x',y')
plot([xw(1,1) xw(1,1)],[y(75) 0],':')
plot([xw(2,1) xw(2,1)],[y(150) 0],':')
plot([xw(3,1) xw(3,1)],[y(250) 0],':')
axis([0 xw(n,1)+.1 -.1 1.1])
text(xw(1,1),-.05,'t_1')
text(xw(2,1),-.05,'t_2')
text(xw(3,1),-.05,'t_3')
text(xw(4,1),-.05,'t_4')
text(xw(5,1),-.05,'t_5')
hold off
