function res = stoch_sim (fun,start,times)

  ## usage:  result = stoch_sim (fun,start,times)
  ##
  ## 

  no_cells=1;
  
  result(1,:)=start;
  len = size(times)(2);
  res = zeros(len,size(start)(2));
  for c=1:no_cells
    c
    for t=2:len
      step=times(1,t)-times(1,t-1);
      dx = feval(fun,result(t-1,:),t);
      result(t,:)=result(t-1,:)+step*dx;
    endfor 
    result;
    res(:,:)+=result(:,:);
  endfor
  res/=no_cells;
endfunction