help with adams.m -> lsode conversion

[Haisam K. Ido]

I'm trying to convert the adams ODE integrator call in matlab 
[t_int,x_int] = adams('xprime',[t_k;t_k1],x_in,1e-10);

In a file called od_main.m (below), to the lsode one below with no success
[t_int,x_int] = lsode( xprime, x_0, [t_k;t_k1] ); 

xprime has this:
function [sys,x0] = xprime(t,x,u,flag)

However, I get the following:

  Initializing variables.
  Loading observation data.
  Iteration number is ...
mloop = 1
error: `flag' undefined near line 11 column 8
error: evaluating expression near line 11, column 8
error: evaluating argument list element number 0
error: evaluating index expression near line 11, column 4
error: evaluating binary operator `==' near line 11, column 13
error: if: error evaluating conditional expression
error: evaluating if command near line 11, column 1
error: called from `xprime' in file
`/disk3/users/idoh/projects/Kalman/kf_labs/od/xprime.m'
error: evaluating argument list element number 0
error: near line 122 of file
`/disk3/users/idoh/projects/Kalman/kf_labs/od/od_main.m'


od_main.m is below:

%
%       OD_MAIN.M
%               Orbit determination using batch least squares.
%

%       Initialize variables.

clear
%close all
format short;

disp('  Initializing variables.')
const;
init;
summary=0;

%       Read in observation data.
         x_k1 = x_k;
         phi_k1 = phi_k;
      else %if nloop~=1,
         x_in = [x_k; reshape(phi_k,225,1)];
%         [t_int,x_int] = adams('xprime',[t_k;t_k1],x_in,1e-10);
        [t_int,x_int] = lsode( xprime, x_0, [t_k;t_k1] );
         [nsize,msize] = size(t_int);
         x_k1 = x_int(nsize,1:15)';
         phi_k1 = reshape(x_int(nsize,16:240)',15,15);
      end

%       Solve for H-tilda and calculated range.

      [h_tilda,rho_calc] = h_calc(x_k1,t_k1,sta_id);

%       Calculate measurement matrix, H, for this iteration.

      h_k1 = h_tilda * phi_k1;

%       Calculate observation error.

      delta_y = rho_k1 - rho_calc;
      summary = [summary
                 delta_y];

%       Accumulate least-square sums.

      sum1 = sum1 + h_k1' * inv(r_obs) * h_k1;
      sum2 = sum2 + h_k1' * inv(r_obs) * delta_y;
      sum_rms = sum_rms + delta_y^2;

%       Update variables for next iteration. End inner loop.

      phi_k = phi_k1;
      x_k = x_k1;
      t_k = t_k1;
   end

%       Compute state error.

   [n_sum,m_sum] = size(sum1);
   [u_svd,s_svd,v_svd] = svd( inv(p_0) + sum1);
   delta_x = zeros(n_sum,1);
   for i_sum=1:n_sum,
      delta_x = delta_x + ( u_svd(:,i_sum)' * sum2 / s_svd(i_sum,i_sum) )
*...
                v_svd(:,i_sum);
   end

%       Check for convergence. (If desired, can keep iterating until
%       rcon is less than some small value.)

   rcon = sqrt(sum_rms/i_size)
   sum_plot = [sum_plot summary];
   summary=[0];

%       Update epoch state.

   x_k = x_0 + delta_x;
   x_0 = x_k;

%       Reset variables for next iteration. End outer loop.

   phi_k = phi_0;
   t_k = t_0;
   sum1 = sum1_init;
   sum2 = sum2_init;
   sum_rms = 0;
end
%end