function [yest,state_vec,t_vec,delta_t] = EKF31(model,param)

Q        =  param.Q;
R        =  param.R;
C        =  param.C;
t_meas   =  param.t_meas;
y_meas   =  param.y_meas;
Endtime  =  param.Endtime;
delta_t  =  param.delta_t;
t_sample = param.t_sample;
simulation_step  =  param.h;

x = param.x0'; % Initial condition 
P = param.P0;  % Initial error covariance

t_est = 0;
j = 1;

% to save estimates
state_vec = zeros(length(t_meas),numel(x));
yest = zeros(numel(t_meas),size(y_meas,2));
t_vec = zeros(length(t_meas),1);

while t_est <= t_meas(end)

    % Measurement update / Correction
    K = P*C'*inv(C*P*C'+R);
    x = x + K*(y_meas(j,:)'-C*x);             % x[k|k]
    P = (eye(4)-K*C)*P;                       % P[k|k]
   
    % Save the state values at t = t_meas
    yest(j,:) = C*x;
    state_vec(j,:) = x;
    t_vec(j) = t_est;
    j = j+1;
    % end
    
    % Time update / Prediction
    [A,x] = feval(model,x,delta_t,t_sample); % Integration step
    P = A*P*A' + Q;                      % P[k+1|k]
    t_est = t_est + t_sample;

end



