clear all;
%% simulation ODE  carModel

t = [0;200];
z0 = [2000; 0; 20];

% solution
[tSol, zSol] = ode45(@ODEcar,t,z0);

figure('Name', 'carSimulation')
subplot(2,2,[1,2]);
plot(tSol, zSol(:,3), 'b-', 'LineWidth', 3);
title('velocity')

subplot(2,2,3);
plot(tSol, zSol(:,2), 'b-', 'LineWidth', 3);
title('position')

subplot(2,2,4);
plot(tSol, zSol(:,1), 'b-','LineWidth', 3);
title('force')

function dzdt = ODEcar(~, z) 
m = 1000;  % mass
k = 50;

F = z(1);
x = z(2);
v = z(3);

% ODE equation 
dxdt = v;
dvdt =  (F - k * v )/ m;

dzdt = [0;dxdt;dvdt];

end