function xhat = EKF(model,init,data)

  M = length(data.y);
  P = init.P0;
  n = numel(init.x0);
  
  xhat = zeros(n,M); xhat(:,1) = init.x0;
  
  for l=1:M-1
    % Measurement update
    xhat(:,l) = FILL_IN_CODE_HERE;

    % Time update
    xhat(:,l+1) = FILL_IN_CODE_HERE;
  end
end