function [y,flag]= romberg(f,a,b,eps1)
%
%  function [y,flag]= romberg(f,a,b,eps1)
%
     format long e;
     h=(b-a)/2;
     n=2;
     c=a+h;
     fa=feval(f,a);
     fb=feval(f,b);
     fc=feval(f,c);
     int=h*(0.5*(fa+fb) +fc)
     coef=4;
     k=1;
     err=1;
     flag=0;
     while err > eps1 & k <= 15
       newval=update(f,a,h,n,int(1))
       for i=1:k
	 temp=int(i);
	 int(i)=newval;
	 newval=(coef(i)*newval - temp)/(coef(i)-1);
       end;
       int = [int newval]
       coef= [coef coef(k)*4 ];
       err=abs(temp-newval)
       k=k+1;
       h=0.5*h;
       n=2*n;
     end;
     y=int(k);
     if err > eps1 
       flag=1;
     end;