How can I improve the speed of my code?

A = 10;
a1 = 1.8;
h1 = 0.4;
b1 = 0.8;
d1 = 2900;
L1 = 0.05;
σ1 = 500;
k = 0.001;
n = 4;

c1[T_] :=
Piecewise[{
{3, 0 < d1*T <= 100},
{2, 100 < d1*T <= 200},
{1,200 < d1*T <= 300},
{0.5, 300 < d1*T},
{0.5, True}}];

TCJS[T_, k_] :=
A/T + c1[T]*d1*T + h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) +
(b1/T)*σ1*Sqrt[T + L1]*
(PDF[NormalDistribution[0, 1], k] -
k*Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}]);

EQ1[T_] :=
(k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/2)*
(PDF[NormalDistribution[0, 1], k] -
k*Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}])*
((T + 2*L1)/(T^2*Sqrt[T + L1])) - A/T^2 + (d1*h1)/2 + c1[T]*d1 == 0;

T1 = T /. Solve[EQ1[T], T, Reals][[1]];

Table[{TCJS[T1, k], T1, k}, {k, 0.001, n, 0.001}]
SortBy[%, First][[1]]


I used Table because I know it is the fastest looping code in Mathematica, but it takes so much more time than matlab.

Is there any way to improve the speed?

• The problem is likely the PDF's, and integrals in your definition of TCJS, because these are being re evaluated every time you evaluate an entry in the table. Using = instead of := when defining TCJS may help. But you'll need to make sure any parameters you don't explicitly list as inputs to the function are set before you define your function. Aug 11 '15 at 15:10
• This code doesn't execute for me - something to do with the Boxes. Could you give us the InputForm of TCJS[T, k] instead? Aug 11 '15 at 15:11
• Take a look at this post on How to copy code so it looks good on this site in order to be able to provide us functional code. Aug 11 '15 at 15:40
• We can reopen your question as soon as you fix the issues pointed out to you, as well as provide information on what you're trying to do. Aug 12 '15 at 1:14
• The code seems fixed now, thanks to m_goldberg. Aug 12 '15 at 14:12

Evaluate the integral once and for all (cf. cdfc):

cdfc[k_] =
Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}];

TCJS[T_, k_] :=
A/T + c1[T]*d1*T +
h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) + (b1/T)*σ1*
Sqrt[T + L1]*(PDF[NormalDistribution[0, 1], k] - k*cdfc[k]);

EQ1[T_] := (k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/
2)*(PDF[NormalDistribution[0, 1], k] -
k*cdfc[k])*((T + 2*L1)/(T^2*Sqrt[T + L1])) -
A/T^2 + (d1*h1)/2 + c1[T]*d1 == 0;

Table[{TCJS[T1, k], T1, k}, {k, 0.001, n, 0.001}]; // AbsoluteTiming


Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

(*  {0.233491, Null}  *)


But we can extend that principal to the other functions as well, provided we block the values of T and k and use the definition of cdfc above:

Block[{T, k},

TCJS[T_, k_] =
A/T + c1[T]*d1*T +
h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) + (b1/T)*σ1*
Sqrt[T + L1]*(PDF[NormalDistribution[0, 1], k] - k*cdfc[k]);

EQ1[T_] = (k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/
2)*(PDF[NormalDistribution[0, 1], k] -
k*cdfc[k])*((T + 2*L1)/(T^2*Sqrt[T + L1])) -
A/T^2 + (d1*h1)/2 + c1[T]*d1 == 0;

]

T1 = T /. Solve[EQ1[T], T, Reals][[1]];

murf = Table[{TCJS[T1, k], T1, k}, {k, 0.001, n, 0.001}]; // AbsoluteTiming


Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

(*  {0.067596, Null}  *)