function [tu,ts,tf,nu,la0]=fopt(name,tu,tf,s0,k,told,tols,mxit,eta)             
%FOPT - Function OPtim. w. terminal constraints, Tf Open
%[tu,ts,tf,nu,la0]=fopt(name,tu,tf,s0,k,told,tols,mxit,eta)
   % Name must be in single quotes; function file 'name' computes 
   % f(s,u)=sdot for flg=1, (Phi,Phis,Phit) for flg=2, [fs,fu] for
   % flg=3; inputs: tu=[t(1:N1),1) u(1:N1,nc)]=init. estim. 
   % u(t); s0(ns,1)=init. state, tf=final time, k=step size param,
   % told=reltol for ode23; stops when norm(dua)<tols or number of
   % iterations > mxit; 0<eta<=1 where d(psi)=-eta*psi is desired
   % change in term. constraints on next iteration; outputs: opt.
   % tu and ts=[t(1:N,1) s(1:N,ns)]), la0=initial adjoint vector.
% BASIC version AEB '84; MATLAB versions Sun H. Hur '90 and
% Fred A. Wiesinger '91, AEB 11/94, 8/97, 2/2/98  
%
if nargin<9, eta=1; end; if nargin<8, mxit=10; end;
[Nu,nc1]=size(tu); nc=nc1-1; it=0; dua=ones(1,nc); 
disp('    Iter.     tf     phi    norm(psi)   dua   dtf');
options=odeset('reltol',told); pmag=1; dtf=1;
while max([norm(dua) pmag abs(dtf)])>tols,
  [t,s]=ode23('fopc_f',[0 tf],s0,options,tu,name); ts=[t s];
  Ns=length(t); ns=length(s0); u1=tu(Nu,[2:nc1]); s1=ts(Ns,[2:ns+1]);
  [Phi,Phis,Phit]=feval(name,u1',s1',tf,2); np=length(Phi);
  phi=Phi(1); psi=Phi([2:np]); pmag=norm(psi)/sqrt(np-1);
  Phid=Phit+Phis*feval(name,u1',s1',tf,1);
  phid=Phid(1); psid=Phid([2:np]);
  yf=[formm(Phis','c'); zeros(np*(np+1)/2,1)]; ny=length(yf);
  [tb,y]=ode23('fopc_b',[tf 0],yf,options,tu,ts,name,ns,np);
  Nb=length(tb); Qg=forms(y(Nb,[np*ns+1:ny])');
  Q=psid*psid'+Qg([2:np],[2:np]); g=psid*phid+Qg([2:np],1);
  nu=-Q\g; la0=formm(y(Nb,[1:np*ns])',ns)*[1; nu];
  du=zeros(Nb,nc); u=du;     
  for i=1:Nb,
     u1=interp1(tu(:,1),tu(:,2:end),tb(i)); 
     s1=interp1(ts(:,1),ts(:,2:end),tb(i));
    % u1=table1(tu,tb(i)); s1=table1(ts,tb(i));
    [fs,fu]=feval(name,u1',s1',tb(i),3);
    la=formm(y(i,[1:np*ns])',ns); HuPhi=la'*fu;
    Huphi=HuPhi(1,:); Hupsi=HuPhi([2:np],:);
    du(i,:)=-k*(Huphi+nu'*Hupsi)-eta*(psi'/Q)*Hupsi;
    u(Nb-i+1,:)=u1+du(i,:);
  end
  dtf=-k*(phid+nu'*psid)-eta*psid'*(Q\psi);
  for j=1:nc, dua(j)=norm(du(:,j))/sqrt(Nb); end;
  if mxit==0, disp([it tf phi pmag dua dtf]); break; end
  disp([it tf phi pmag dua dtf]);
  tf=tf+dtf; t1=zeros(Nb,1); for i=1:Nb, t1(i)=tb(Nb-i+1); end;
  tu=[t1*tf/t1(Nb) u]; Nu=Nb; tu(Nb,1)=tf*(1+1e-8);
  if it>=mxit, break, end; it=it+1;
  end
   
   