need help to calculate a dependent variable outside of function using ode45

Below is my function and the code to call that function using ode45 for a first order non-homogenius equation. THe code works fine , and the answer is also correct. But i just want to calculate kf outside of the program, to decrease the runtime

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = FunctionOfNacheilzone(N,q)                   % Function to calculate stress in lagging zone
 
Radius_R = 21;                                            % Radius of  idle Roll 
n = 1;
%___________________________________________________________________________
 
Wire_dH = 8;                                              % Diameter of hot rolled copper                                                             
Wire_d0 = 2;                                              % Reduced diameter for roll drwaing operation
meu = 0.1;                                                % Friction acting in between wire and rolls
C = 522;                                                  % Coefficient of Ludwik equation                               
M = 0.368;                                                % Coefficiant of Ludwik equation
%___________________________________________________________________________
                                   
h_1 = [1.7900 ;  1.7000;    1.6000;    1.5000;    1.4000]; % Final height after reduction                     
h_0 = [2; 2; 2; 2; 2 ];                                   % Initi height of copper wire                   
b_1 = [2.0400;   2.0900;  2.1400;  2.2100;    2.2700];    % Final width after roll drawing 
b_0 = [2; 2; 2; 2; 2 ];                                   % Initial width of copper wire
Del_h = abs(h_1(n,1) - h_0(n,1));                         % Relative height reduction
l_d = sqrt(Radius_R.*Del_h);                              % Bite angle 
b = ((b_1(n,1)-b_0(n,1))./l_d)*(l_d-N) +b_0(n,1);               % Width as a function of position shift in XY plane
h = h_1(n,1)+(N^2/Radius_R);                              % Height as a function of position shift in XY plane  
pi_b = log(b./b_0(n,1));                                  % Degree of deformation along with width
pi_h = log (h./h_0(n,1));                                 % Degree of deformation along with height
pi_l = (-pi_b - pi_h);                                    % Degree of deformation along with length
pi_eq  = sqrt((2/3).*(((pi_l.^2)+(pi_b.^2)+(pi_h.^2))));  % Equivalent degree of deformation
pi_0 = abs(log(Wire_d0^2/Wire_dH^2));                     % Initial degree of deformation
pi_1 = pi_eq + pi_0;                                      % Increase of degree of deformation with increase of span
kf = C.*(pi_1).^(M);                                      % Kf values at every point by Ludwik equation
%_________________________________________________________________________
% the ordinary differential equation of first order to determine stress in
% the lagging zone 
dqdN = -(((2./(h_1+N^2/Radius_R)).*((tan(N/Radius_R)-tan((N/Radius_R)-atan(meu))))).*q)...
    -((2./(h_1+N^2/Radius_R)).*kf.*tan((N/Radius_R)-atan(meu)));
f = dqdN; 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 1;                                             
Radius_R = 21;                                     % Radius of idle roll              
h_1 = h1(1:end,:);                                 % Height after reduction
h_0 = h0(1:end,:);                                 % Initial height of copper wire
Delta_h = abs(h_1(n,1) - h_0(n,1));                % Reltive height reduction
In_alf = 0;                                        % Final condition of bite angle for lagging zone
Fi_alf = sqrt(Radius_R.*Delta_h);                  % Initial condition of bite angle
Nspan = [Fi_alf:-0.0001:In_alf];                   % Span of bite angle
 
%____________________________________________________

q0 = [0;0;0;0;0];                                % Intital condition of stress in the lagging zone for ODE
[Nout, q] = ode45(@FunctionOfNacheilzone,Nspan,q0); % Calling the ode function using MATLAB built-in function

Answer 1

[Walter Roberson]

kf = 522*(((2*(log(51/50 - N/105) + log(N^2/42 + 179/200))^2)/3 + (2*log(N^2/42 + 179/200)^2)/3 + (2*log(51/50 - N/105)^2)/3)^(1/2) + 6243314768165359/2251799813685248)^(46/125)

That needs N to calculate, so you cannot pre-calculate it.

A few parts of it can be pre-calculated, but not so much.

For whatever can be pre-calculated, see http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[md mayeen uddin]

sorry for a slight mistake , N = Nspan(i have corrected the code). so that is already a defined range.Can it now be solvable outside of function?

[Walter Roberson]

No, that does not help.

syms N
syms q [5 1]
Radius_R = 21;                                            % Radius of  idle Roll 
n = 1;
%___________________________________________________________________________
 
Wire_dH = 8;                                              % Diameter of hot rolled copper                                                             
Wire_d0 = 2;                                              % Reduced diameter for roll drwaing operation
meu = 0.1;                                                % Friction acting in between wire and rolls
C = 522;                                                  % Coefficient of Ludwik equation                               
M = 0.368;                                                % Coefficiant of Ludwik equation
%___________________________________________________________________________
                                   
h_1 = [1.7900 ;  1.7000;    1.6000;    1.5000;    1.4000]; % Final height after reduction                     
h_0 = [2; 2; 2; 2; 2 ];                                   % Initi height of copper wire                   
b_1 = [2.0400;   2.0900;  2.1400;  2.2100;    2.2700];    % Final width after roll drawing 
b_0 = [2; 2; 2; 2; 2 ];                                   % Initial width of copper wire
Del_h = abs(h_1(n,1) - h_0(n,1));                         % Relative height reduction
l_d = sqrt(Radius_R.*Del_h);                              % Bite angle 
b = ((b_1(n,1)-b_0(n,1))./l_d)*(l_d-N) +b_0(n,1);               % Width as a function of position shift in XY plane
h = h_1(n,1)+(N^2/Radius_R);                              % Height as a function of position shift in XY plane  
pi_b = log(b./b_0(n,1));                                  % Degree of deformation along with width
pi_h = log (h./h_0(n,1));                                 % Degree of deformation along with height
pi_l = (-pi_b - pi_h);                                    % Degree of deformation along with length
pi_eq  = sqrt((2/3).*(((pi_l.^2)+(pi_b.^2)+(pi_h.^2))));  % Equivalent degree of deformation
pi_0 = abs(log(Wire_d0^2/Wire_dH^2));                     % Initial degree of deformation
pi_1 = pi_eq + pi_0;                                      % Increase of degree of deformation with increase of span
kf = C.*(pi_1).^(M);                                      % Kf values at every point by Ludwik equation

simplify(kf)

ans = 

%_________________________________________________________________________
% the ordinary differential equation of first order to determine stress in
% the lagging zone 
dqdN = -(((2./(h_1+N^2/Radius_R)).*((tan(N/Radius_R)-tan((N/Radius_R)-atan(meu))))).*q)...
    -((2./(h_1+N^2/Radius_R)).*kf.*tan((N/Radius_R)-atan(meu)));
f = dqdN; 

f

f = 

Your kf depends upon N, which is the current "time" (or at least the independent variable.) Pre-computing would only work if kf was independent of the inputs N and q, or if it depended on N and q in "simple" ways so that you could break out part of the computation . But kf does not depend upon N and q in simple ways: that entire sqrt() expression needs to be evaluated according to the current value of the independent variable.

Please also notice the in the above output of f, which has . If 51/50 - N/105 is negative then the log() will be complex. Can we prove that N/105 will be less than 51/50 ?

If we make guesses that h1 is the same as h_1 and h0 is the same as h_0 then

Radius_R = 21;                                     % Radius of idle roll              
Delta_h = abs(h_1(n,1) - h_0(n,1));                % Reltive height reduction
In_alf = 0;                                        % Final condition of bite angle for lagging zone
Fi_alf = sqrt(Radius_R.*Delta_h);                  % Initial condition of bite angle
Nspan = [Fi_alf:-0.0001:In_alf];                   % Span of bite angle
max(Nspan)/105

ans = 0.0200

In this case the answer is Yes, N/105 is comfortably less than 51/50 ... provided that our guess was right about h1 vs h_1 and h0 vs h_0