Need help about simulation of vehicle dynamics

Hi all,

I recently have been trying to reconstruct the steps described in the paper 'Driving with Tentacles - Integral Structures for Sensing and Motion'. I want to recreate its simulation of vehicle dynamics in MATLAB and I really want to know if I'm doing it correctly or not. I also would like a detailed explanation about the simulation method described in the paper. Below is my code in 'vehicleSim.m', a thousand thanks in advance!

function [  ] = vehicleSim( vi )
% initial conditions
m = 2900; % (kg) car weight
Jz = 5561; % (kgm.^2) moment of inertial of the car
caf = 80000; % (N/rad) cornering stiffness of front axis
car = 110000; % (N/rad) cornering stiffness of rear axis
lf = 1.425; % (m) distance from front axis to center of gravity
lr = 1.425; % (m) distance from rear axis to center of gravity
dt = 0.001; % (s) integration step
ecp = 0.2; % (m) distance from wind pressure point to center of gravity
Fxr = 0; % (N) longitudal force on rear axis
Fxf = 0; % (N) longitudal force on front axis
Fax = 0; % longitudal wind force
Fay = 0; % lateral wind force
beta = 0; % sideslip angle
phiDot = 0; % yaw angle rate
v = vi; % (m/s) initial speed
df = 2;
beta_set = [];
phiDot_set = [];
v_set = [];
t_set = [];
x_set = [];
y_set = [];
r_set = [];
x = 0;
y = 0;
phi = 0; % yaw angle
r_sum = 0;
% r_target = 4.244131815783875;
% v = 0.25;
% for deltaF = -0.5*pi:df:0.5*pi
for deltaF = -80:df:80
  deltaF = deltaF/180*pi
  for i = 0:dt:50
      % state transit
      Fyr = car * atan((lr * phiDot - v * sin(beta))/v/cos(beta));
      Fyf = caf * (deltaF - atan((lf * phiDot + v * sin(beta))/v/cos(beta)));

        betaDot = -phiDot - sin(beta) * (Fxr - Fax + Fxf * cos(deltaF) - Fyf * sin(deltaF)) /m/v - cos(beta) * (-Fyr - Fay - Fxf * sin(deltaF) - Fyf * cos(deltaF))/m/v;
        phiDDot = 1/Jz*(Fyf * cos(deltaF) + Fxf * sin(deltaF))*lf - Fyr * lr + Fay * ecp;
        vDot = cos(beta) * (Fxr - Fax + Fxf*cos(deltaF) - Fyf*sin(deltaF))/m - sin(beta) * (-Fyr - Fay - Fxf*sin(deltaF) - Fyf*cos(deltaF))/m;

        if betaDot < 0.01 || phiDDot < 0.01 || vDot < 0.01
            i,deltaF
            break;
        end

        % get current state
        beta = beta + betaDot*dt;
        phiDot = phiDot + phiDDot*dt;
%         phiDot = v/r_i;
        phi = phi + phiDot*dt;
        v = v + vDot*dt;
        x = x + v*cos(beta + phi)*dt;
        y = y + v*sin(beta + phi)*dt;        
%         deltaF = atan(phiDot*(lf+lr)/v);     
        beta_set = [beta_set beta];
        phiDot_set = [phiDot_set phiDot];
        v_set = [v_set v];
        t_set = [t_set i]; 
        x_set = [x_set x];
        y_set = [y_set y];
        r_set = [r_set v/phiDot];
        r_sum = r_sum+v/phiDot;

    end
    % plot
    figure

    subplot(3,1,1);
    plot(t_set,v_set,'.')
    str_vdeltaF = num2str(deltaF/pi*180);
    str_r = num2str(r_sum/size(r_set,2));
    acc_str = num2str(v*v/phiDot);
    v_final_str = num2str(v);
    beta_final = num2str(beta);
    title(['deltaF: ',str_vdeltaF,' r: ',str_r,' acc: ',acc_str,' vFinal: ',v_final_str,' betaFinal: ',beta_final])   

    subplot(3,1,2);
    plot(t_set,r_set,'.')
    title(['radius'])

    subplot(3,1,3);
    plot(x_set,y_set,'.');
    title(['coords'])
    axis equal

    % reinitialize
    beta = 0;
    phiDot = 0;
    v = vi; % (m/s)
    x = 0;
    y = 0;
    phi = 0;
    r_sum = 0;
    beta_set = [];
    phiDot_set = [];
    v_set = [];
    t_set = [];
    x_set = [];
    y_set = [];
    r_set = [];
end

end