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. – N.J.Evans Aug 11 '15 at 15:10 4   This code doesn't execute for me - something to do with the Boxes. Could you give us the InputForm of TCJS[T, k] instead? – Patrick Stevens 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. – MarcoB 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. – J. M.♦ Aug 12 '15 at 1:14      The code seems fixed now, thanks to m_goldberg. – Michael E2 Aug 12 '15 at 14:12 ================= 1 Answer 1 ================= 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} *)