function [tu,ts,la0]=fop0(name,tu,tf,s0,k,told,tols,mxit)             
%FOP0 - Function OPtim. w. 0 terminal constraints
%[tu,ts,la0]=fop0(name,tu,tf,s0,k,told,tols,mxit)
   % Name must be in single quotes; function file 'name' computes 
   % sdot=f(s,u) for flg=1, (phi,phis) for flg=2, fs for flg=3,
   % fu for flg=4; inputs: tu=[t(N1),1) u(N1,nc)]=initial estimate
   % of u(t); s0(ns,1)=initial state, tf=final time, k=step size
   % parameter, told=reltol for ode23; stops when norm(dua)<
   % tols or no. iterations > mxit; outputs: optimal tu and ts 
   % where ts=[t(N,1) s(N,ns)]), la0=initial adjoint vector;
% BASIC version AEB '84; MATLAB versions Sun H. Hur '90, Fred
% A. Wiesinger '91, Paul M. Montgomery '95, AEB 11/94, 2/2/98  
%
% modifed 3/21/0b by R. Throne, changed table1 to interp1
%
if nargin<8, mxit=10; end; it=0; dua=1; 
disp('     Iter.       phi        dua'); 
options=odeset('reltol',told); [dum,nc1]=size(tu); nc=nc1-1;
dum=zeros(nc,1);
while norm(dua)>tols, 
   [t,s]=ode23('fop0_f',[0 tf],s0,options,tu,name); ts=[t s];
   Ns=length(t); [phi,phis]=feval(name,dum,s(Ns,:)',tf,2); 
   [tb,la]=ode23('fop0_b',[tf 0],phis',options,tu,ts,name);
   N=length(tb); la0=la(N,:); u=zeros(N,nc); du=u; 
   for i=1:N, 
     s1=interp1(ts(:,1), ts(:,2:end),tb(i)); u1=interp1(tu(:,1),tu(:,2:end),tb(i));
     [fs,fu]=feval(name,u1,s1,tb(i),3); Hu=la(i,:)*fu;
     du(i,:)=-k*Hu; 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 dua]); break; end
   disp([it phi 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
