function [tu,ts,nu,la0]=fopc(name,tu,tf,s0,k,told,tols,mxit,eta)             
%FOPC - Function OPtim. w. terminal Constraints
%[tu,ts,nu,la0]=fopc(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) for flg=2, [fs,fu] for flg=3,
   % inputs: global variables name, tu, ts, ns, np, where np-1 = no.
   % terminal constraints, tu=[t(N1),1) u(N1,nc)]= init. estim. of
   % 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(N,1) s(N,ns)]), nu and la0=initial adjoint vector.
% BASIC version AEB '84; MATLAB versions Sun H. Hur '90 and
% Fred A. Wiesinger '91, AEB 11/94, 2/2/98  
%
if nargin<9, eta=1; end; if nargin<8, mxit=10; end;
it=0; dua=1; pmag=1;         
disp('     Iter.       phi        norm(psi)     dua');
options=odeset('reltol',told); [dum,nc1]=size(tu); nc=nc1-1; 
while max([norm(dua) pmag])>tols, 
  [t,s]=ode23('fopc_f',[0 tf],s0,options,tu,name); ts=[t s];
  Ns=length(t); ns=length(s0); dum=zeros(nc,1); 
  [Phi,Phis]=feval(name,dum,s(Ns,:)',tf,2); np=length(Phi);
  phi=Phi(1); psi=Phi([2:np]); pmag=norm(psi)/sqrt(np-1);
  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);
  N=length(tb); Qg=forms(y(N,[np*ns+1:ny])'); g=Qg([2:np],1);
  Q=Qg([2:np],[2:np]); nu=-Q\g; 
  la0=formm(y(N,[1:np*ns])',ns)*[1; nu]; du=zeros(N,nc); u=du;
  for i=1:N, 
    u1=interp1(tu(:,1),tu(:,2:end),tb(i)); s1=interp1(ts(:,1),ts(:,2:end),tb(i));
    % s1=table1(ts,tb(i)); u1=table1(tu,tb(i)); modified by rdt 3/29/06
    [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(N-i+1,:)=u1+du(i,:);
  end
  for j=1:nc, dua(j)=norm(du(:,j))/sqrt(N); end;
  if mxit==0, disp([it phi pmag dua]); break; end
  disp([it phi pmag dua]);
  t1=zeros(N,1); for i=1:N, t1(i)=tb(N-i+1); end; tu=[t1 u];
  if it>=mxit, break, end; it=it+1;
end
   
   