I am trying to build a PID controller for a system with Matlab code . If i have to design it using the Eulers method, how to proceed

Question

I have a third order maglev model for which I want to design a PID controller. I wish to build the controller using Euler's loop considering the system dynamics. The purpose of building the controller is to stabilize the ball position. The states are taken as x1,x2,x3 as position, velocity and current. Please suggest how to build the matlab code for the same.

I have tried the following code.

% Constants
g       = 9.81;
k1      = 2.48315625e-5;
R       = 2.52;
L       = 0.277;
a       = 143.48;
b       = -2.8;
m       = 0.02;
l       = 0;
u1      = 0; % Step size
h       = 0.001; % Time step size
 x1      = 0.0125; % Ball Position(meter) initial value
 x2      = 0; % Velocity of the Ball
 x3      = 1.111; % Initial current

x1d      = 0.009;   % Desired position in meters
x3d     = 0.8;      % Desired current
e1      = x1 - x1d; % Position error
e2      = 0;        % Velocity error
e3      = x3 - x3d; % Current error

% Initial values for simulation

 x1void   = 0.0125; 
x2void  = 0;
x3void  = 1.111;
e1void  = e1; 
e2void  = 0;
e3void  = e3;
t       = 0;

% PID controller parameters
Kp      = 4; % Proportional gain
Ki      = 2; % Integral gain
Kd      = .2; % Derivative gain

% Initialize PID components
integral_error = 0; 
previous_error = e1; 

% e1      = e1 / a; % Convert voltage to meters (a=143.48)
 u = 0; % Control input initialization

% Simulation loop
 while t <= 10
% Error derivatives for PID controller
e1_dot = e2; % Derivative of position error
e2_dot = g - (k1 * (0.8 + e3) * (0.8 + e3)) / (m * (0.009 + e1) * (0.009 + e1)); % Velocity   equation
e3_dot = -(R / L) * (0.8 + e3) + (2 * k1 * e2 * (0.8 + e3)) / (L * (0.009 + e1) * (0.009 + e1)) + (1 / L) * u; % Current equation

% PID controller calculation
error = e1;  % Current position error
integral_error = integral_error + error * h;  % Integral of error
derivative_error = (error - previous_error) / h;  % Derivative of error
u = Kp * error + Ki * integral_error + Kd * derivative_error;  % PID control law

% Update state variables
e1 = e1 + e1_dot * h;
e2 = e2 + e2_dot * h;
e3 = e3 + e3_dot * h;

% Update previous error
previous_error = error;

% Time update
t = t + h; 
l = [l; t];
u1 = [u1, u];

% Ball dynamics
x = e1 + x1d; % Ball position
x2 = e2; % Velocity of the ball
x3 = e3 + x3d; % Current
x1 = a * x + b; % Sensor output voltage

% Store variables for plotting
e1void = [e1void; e1];
e2void = [e2void; e2];
e3void = [e3void; e3];
x1void = [x1void; x1];
x2void = [x2void; x2];
x3void = [x3void; x3];

end

% Plot results
figure(1);
plot(l, x1void);
grid on;
% axis([0 10 0.0087 0.0125]);
title('Ball Position (meter)');
xlabel('Time (sec)');
ylabel('Ball Position (meter)');

figure(2);
plot(l, x2void);
grid on;
% axis([0 10 -0.005 0.001]);
title('Ball Velocity');
xlabel('Time (sec)');
ylabel('Ball Velocity (m/s)');

figure(3);
plot(l, x3void);
grid on;
% axis([0 10 0.75 1.15]);
title('Current');
xlabel('Time (sec)');
ylabel('Current (A)');

figure(4);
plot(l, u1);
grid on;
axis([0 10 min(u1)-1 max(u1)+1]);
title('Control Input');
xlabel('Time (sec)');
ylabel('Control Input (V)');