How to merge external feedback control (function) into the original S-function (system model)?

This is my original level 1 s-function:

   function [sys,x0,str,ts] = LongTyre(t,x,u,flag)

% definition of constants
m      = 440.00;         % [kg]
I      =  40.00;         % [kg*m2]
R      =   0.32;         % [m]
Fn0    =   9.81*m;       % [N] 
v0     =  40.0;          % [m/s] 
vr0    =   0.0;          % [m/s]
w0     =  1/R*(v0-vr0);  % [-/s]

% LuGre friction model parameters
alpha  =    2.0;         % [-]
sigma0 =  181.54;        % [1/m]
sigma1 =    1.00;        % [1/m]
sigma2 =    0.0018;      % [s/m] 
muC    =    0.8;         % [-]
muS    =    1.55;        % [-] 
vS     =    6.57;        % [m/s]
vd     =    0.001;       % [m/s]
L      =    0.2;         % [m]
kappa  =    1.6667*L;    % [-]
theta  =    0.2;         % [-]

% Initial bristle deflection value 
h0    = theta*(muC + (muS-muC)*exp(-abs(vr0/vS)^alpha));
k00   = sigma0*Fn0;
g00   = abs(vr0)/h0/Fn0;
g10   = kappa*R*abs(w0)/k00;
gt0   = g00+g10;
if vr0==0;
    z0 = 0;
else
    z0 = vr0/gt0/k00;
end

switch flag,

    %%%%%%%%%%%%%%%%%%
    % Initialization %
    %%%%%%%%%%%%%%%%%%
    case 0     
      [sys,x0,str,ts] = mdlInitializeSizes(z0,v0,w0);                                

    %%%%%%%%%%%%%%%
    % Derivatives %
    %%%%%%%%%%%%%%%
    case 1,
      sys = mdlDerivatives(t,x,u,m,R,I,alpha,sigma0,sigma1,sigma2,muC,muS,vS,vd,L,kappa,theta,Fn0);

    %%%%%%%%%%%
    % Outputs %
    %%%%%%%%%%%
    case 3,
      sys = mdlOutputs(t,x,u,m,R,I,alpha,sigma0,sigma1,sigma2,muC,muS,vS,vd,L,kappa,theta,Fn0);

    %%%%%%%%%%%%%
    % Terminate %
    %%%%%%%%%%%%%
    case 9                                                
      sys = []; % do nothing

end

%end simom

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes(z0,v0,w0)

sizes = simsizes;

sizes.NumContStates  = 3;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 10;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 2;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);
x0  = [z0,v0,w0]; %zeros(sizes.NumContStates,1); 
str = [];
ts  = [0 0];

% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys = mdlDerivatives(t,x,u,m,R,I,alpha,sigma0,sigma1,sigma2,muC,muS,vS,vd,L,kappa,theta,Fn0)

% rename state variables
z  = x(1);
v  = x(2);
w  = x(3);

% rename inputs
T  = u(1);
Fn = u(2);

% equations of motion (v_rel=v-R*w and Fw=-alpha*v_rel)
vr   = (v-R*w);
g    = theta*(muC+(muS-muC)*exp(-abs(vr/vS)^alpha));
k0   = sigma0*(Fn+Fn0);
g0   = abs(vr)/g/(Fn+Fn0);
g1   = kappa*R*abs(w)/k0;
gt   = g0+g1;
dzdt = vr-k0*gt*z;   
F    = -(sigma0*z + sigma1*dzdt + sigma2*vr)*(Fn0+Fn);

dvdt = 1/m*F;
dwdt = 1/I*(-F*R+T);

% rename state derivatives of variables
sys(1) = dzdt;
sys(2) = dvdt;
sys(3) = dwdt;

% end mdlDerivatives
%
%=============================================================================
% mdlOutputs
% Return the output vector for the S-function
%=============================================================================
%
function sys = mdlOutputs(t,x,u,m,R,I,alpha,sigma0,sigma1,sigma2,muC,muS,vS,vd,L,kappa,theta,Fn0)

% rename state variables
z  = x(1);
v  = x(2);
w  = x(3);

% rename inputs
T  = u(1);
Fn = u(2);

% output definition
if and(v==0,w==0);
    s = 0;
elseif abs(R*w)>abs(v);
    s = 1-v/(R*w);
else;
    s = R*w/v-1;
end;

% friction force
vr   = (v-R*w);
g    = theta*(muC + (muS-muC)*exp(-abs(vr/vS)^alpha));
k0   = sigma0*(Fn+Fn0);
g0   = abs(vr)/g/(Fn+Fn0);
g1   = kappa*R*abs(w)/k0;
gt   = g0+g1;
zdot = vr-k0*gt*z;
F    = -(sigma0*z + sigma1*zdot + sigma2*vr)*(Fn0+Fn);

% output
sys = [ z v w F T vr gt s zdot k0 ]; 

% end mdlOutputs

Using an external control law like this:

function T=abs_u(z,vr,Fn)

% Model Parameters
m      = 440.00;         % [kg]
I      =  40.00;         % [kg*m2]
R      =   0.32;         % [m]
Fn0    =   9.81*m;       % [N] 
mr=(1/m+R^2/I)^-1;       % [kg]

% LuGre friction model parameters
sigma0 = 181.54;         % [1/m]
muS    = 1.55;           % [-]
muSh   = muS+1;

% Control Optional Parameters
a      = 0.5;            % [-] % Assumed to be 0.5;
zs=0.084;
vrs=2;

k0     = sigma0*(Fn+Fn0);

f      = z+vr^2/(2*muSh);
fs     = zs+vrs^2/(2*muSh);

T=I/mr/R*(2*a*(f-fs)-k0*zs);

Takes minutes to solve this (above 60mins!!!) How can I merge the control law in the original s-function? Is it even possible?

Note: The u(1)=T is control law while u(2)=Fn is external input.