Electro-hydraulic actuator system simulating based on ode object generates error

I want to simulate the electro-hydraulic actuator (EHA) system using MATLAB's ode object. The schematic diagram of the EHA is shown as follows.

However, the code does not generate effective result. I think the error is induced by the stiffness of the EHA system, what should I do to obtain the effective result?

% EHA system parameters initialization
EHAParams = struct;
EHAParams.A1 = 1.96e-3; % Hydraulic cylinder piston side area
EHAParams.A2 = 1.65e-4; % Hydraulic cylinder rod side area
EHAParams.L = 0.5; % Cylinder stroke
EHAParams.V10 = EHAParams.A1 * EHAParams.L / 2; % Initial Volume V10 [mid position]
EHAParams.V20 = EHAParams.A2 * EHAParams.L / 2; % Initial Volume V20 [mid position]

EHAParams.beta = 1700e6; % Effective bulk modulus
EHAParams.b = 100; % Friction coefficient
EHAParams.k = 20; % Stiffness coefficient
EHAParams.m = 50; % Load mass
EHAParams.Pt = 0; % Tank pressure
EHAParams.Ps = 21e6; % Supply pressure
EHAParams.FL = 100; % Constant load

EHAParams.rho = 850; % Fluid density
EHAParams.Cd = 0.62; % Discharge coefficient
EHAParams.Wv = pi*1e-3/4; % Servo valve area gradient
EHAParams.ku = 2e-4; % Spool displacement gain

% PID controller parameters
 EHAParams.PID.Kp = 200;
 EHAParams.PID.Ki = 50;
 EHAParams.PID.Kd = 1; 
 % reference position and velocity cmd for PID
 EHAParams.refPos = @(t) 0.1 + 0.05*sin(pi*t);
 EHAParams.refVel = @(t) 0.05*pi*cos(pi*t);

% Simulate the EHA system using ode object
X0 = [0;0;1e7;1e7;0]; % initial state [pos;vel;P1;P2;error_int]
startTime = 0;
endTime = 10;

F = ode(ODEFcn=@(t,x,p) stateTransitionEHA(t,x,p),InitialTime=startTime,InitialValue=X0,...
        Parameters=EHAParams,Solver="ode23t");
F.SolverOptions.OutputFcn = @odeplot;
F.SolverOptions.OutputSelection = [3 4];

sol = solve(F,startTime,endTime);

Error using odePlotImpl (line 61)
Error updating the ODEPLOT window. Solution data may have been corrupted. Argument Y cannot be complex.

Error in odeplot (line 40)
status = odePlotImpl(fun,flag,"odeplot",SetAxis=setAxis);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ode23t (line 758)
stop = feval(outputFcn,tout_new,yout_new(outputs,:),'',outputArgs{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.ode.internal.LegacySolver/runSolver (line 672)
[varargout{1:nout}] = obj.SolverFcn(obj.ODEFcn, tspan, ...
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in matlab.ode.internal.LegacySolver/solve (line 141)
[t,y] = obj.runSolver(obj.InitialTime,t2);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ode/solveBetween (line 1009)
[tR,yR] = solver.solve(t0,tmax,pps);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ode/solve (line 388)
sol = obj.solveBetween(t1,t2,refine);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function dX = stateTransitionEHA(t,X,Params)
    % State transition model of EHA 
    % Input arguments: t->current time
    %                  X->current state                 
    %                  Params -> struct contains EHA parameters
    % Output arguments: dX -> state differential

     H = @(u) max(0,sign(u)); % Heaviside function H
     
     dX = zeros(5,1); % Preallocation state derivative

     dX(1) = X(2);
     dX(2) = (-Params.k*X(1) - Params.b*X(2) + Params.A1*X(3) - Params.A2*X(4) - ...
               Params.FL)/Params.m; % 
     u = Params.PID.Kp*(Params.refPos(t) - X(1)) + Params.PID.Kd*(Params.refVel(t) - X(2)) + ...
         Params.PID.Ki*X(5); % current control input
     % ------------------------------------------------------- %
     kq = Params.Cd * Params.Wv * sqrt(2/Params.rho); 
     g1 = H(u)*sqrt(Params.Ps - X(3)) + H(-u)*sqrt(X(3));
     g2 = H(u)*sqrt(X(4)) + H(-u)*sqrt(Params.Ps - X(4));
     h1 = 1 / (Params.V10 + Params.A1 * X(1));
     h2 = 1 / (Params.V20 - Params.A2 * X(1));

     dX(3) = Params.beta * h1 * (kq * Params.ku * g1 * u - Params.A1 * X(2) - ...
                                 2e-6 * (X(3) - X(4)));
     dX(4) = Params.beta * h2 * (-kq * Params.ku * g2 * u + Params.A2 * X(2) + ...
                                 2e-6 * (X(3) - X(4)));
     dX(5) = Params.refPos(t) - X(1); % tracking error (ref - pos)
end