We are facing a problem in the following code where we intend to solve set of differential equations using Range Kutta method of order 4.

function w = rk4_systems(a,b,N,init) format longG init = [945.027,0,0,0,0,0,0,0,0,0,0,0,0] global E1 E3 E10 E11 E17 E19 E20 E21 E22 K01 K03 K010 K011 K017 K019 K020 K021 K022 R E1 = 52.58; E3 = 65.33; E10 = 64.17; E11 = 56.9; E17 = 50.73; E19 = 34.56; E20 = 35.64; E21 = 57.97; E22 = 29.76; K01 = 6.565E11; K03 = 7.284E12; K010 = 7.386E12; K011 = 2.424E11; K017 = 2.075E11; K019 = 8.385E12; K020 = 9.74E11; K021 = 6.4E17; K022 = 1.51E12; R = 1.98E-3; D = 5; L = 9000; V = pi*(D^2)*L/4; m = size(init,1); if m == 1 init = init'; end h = (b-a)/N; %the step size x(1) = a; w(:,1) = init; %initial conditions

for i = 1:N k1 = h*f(x(i), w(:,i)); k2 = h*f(x(i)+h/2, w(:,i)+0.5*k1); k3 = h*f(x(i)+h/2, w(:,i)+0.5*k2); k4 = h*f(x(i)+h, w(:,i)+k3); w(:,i+1) = w(:,i) + (k1 + 2*k2 + 2*k3 + k4)/6; x(i+1) = a + i*h;

   T1 = -0.961*x(i)^2 + 28.667*x(i) + 862.16;
    T2= -0.961*(x(i-1))^2 + 28.667*(x(i-1)) + 862.16;
     T = (T1+T2)/2;
     T(:,1)=T;
    K1 = K01*(exp(-E1/(R*T)));
    K3 = K03*(exp(-E3/(R*T)));
    K10 = K010*(exp(-E10/(R*T)));
    K11 = K011*(exp(-E11/(R*T)));
    K17 = K017*(exp(-E17/(R*T)));
    K19 = K019*(exp(-E19/(R*T)));
    K20 = K020*(exp(-E20/(R*T)));
    K21 = K021*(exp(-E21/(R*T)));
    K22 = K022*(exp(-E22/(R*T)));

    M1 = [1.79E-7 0.685 -0.59]
M2 = [5.25E-7 0.59006 105.67]
M3 = [2.0789E-6 0.4163 352.7]
M4 = [7.391E-7 0.5423 263.73]
M5 = [6.97E-7 0.5462 305.25]
M6 = [6.0259E-7 0.5309 199.64]
M7 = [1.7514E-7 0.70737 157.14]
M8 = [3.13E-8 0.9679 7.9]
M9 = [8.7268E-7 0.49397 323.79]
M10 = [8.3436E-7 0.49713 365.86]
M11 = [6.3863E-7 0.5254 295.1]
M12 = [4.22E-7 0.58154 239.21]
Mhept = [6.77E-8 0.8367 36.7]
Mch = [6.672E-8 0.82837 85.752]
Msteam = [1.7096E-8 1.146 0]

Mu1 = (M1(1,1)*(T)^(M1(1,2)))/(1+(M1(1,3)/T)) Mu2 = (M2(1,1)*(T)^(M2(1,2)))/(1+(M2(1,3)/T)) Mu3 = (M3(1,1)*(T)^(M3(1,2)))/(1+(M3(1,3)/T)) Mu4 = (M4(1,1)*(T)^(M4(1,2)))/(1+(M4(1,3)/T)) Mu5 = (M5(1,1)*(T)^(M5(1,2)))/(1+(M5(1,3)/T)) Mu6 = (M6(1,1)*(T)^(M6(1,2)))/(1+(M6(1,3)/T)) Mu7 = (M7(1,1)*(T)^(M7(1,2)))/(1+(M7(1,3)/T)) Mu8 = (M8(1,1)*(T)^(M8(1,2)))/(1+(M8(1,3)/T)) Mu9 = (M9(1,1)*(T)^(M9(1,2)))/(1+(M19(1,3)/T)) Mupi = (0.7*Mu8) + (0.2*Mu9) + (0.1*Mu7) Mu10 = (M10(1,1)*(T)^(M10(1,2)))/(1+(M10(1,3)/T)) Mu11 = (M11(1,1)*(T)^(M11(1,2)))/(1+(M11(1,3)/T)) Mu12 = (M12(1,1)*(T)^(M12(1,2)))/(1+(M12(1,3)/T)) Muhept = (Mhept(1,1)*(T)^(Mhept(1,2)))/(1+(Mhept(1,3)/T)) Much = (Mch(1,1)*(T)^(Mch(1,2)))/(1+(Mch(1,3)/T)) Mnaphtha = (Muhept+Much)/2 end

[x' T' w'] function dC = f(T,C) dC = [-(K1*C0);(2*K19*C3*C6)+(2*K20*C6*C4)+(2*K21*C6*C5)+(2*K22*C6^2)+(K1*0.58*C0);(K3*C4)-(3*K11*C4^2)+(0.68*K1*C0);-(K3*C4^2)-(K19*C3*C6)+(0.88*K1*C0);-(K3*C4)-(2*K10*C4^2)-(2*K11*C4^2)-(K20*C6*C3)+(0.6*K1*C0);-(K1*0.2*C0)-(K17*C5)-(K21*C6*C5);(0.07*K1*C0)-(K19*C3*C6)-(K20*C6*C4)-(K21*C6*C5)-(2*K22*C6^2);(K11*0.14*C4^2)+(0.19*K17*C5);(K19*C6*C3)+(0.3*K11*C4^2)+(0.41*K17*C5);(K20*C5*C4)+(K11*C4^2)+(0.41*K17*C5);(K11*C4^2)+(0.41*K17*C5);(K22*C5^2);(K21*C6*C5)];

Answer 1

[Jan]

The code is not formatted properly - unfortunately. Please care for posting coede in a form, which allows to read it, when you expect others to do so.

But here the problem is easy: You cannot create functions inside a script or in the command window. An M-file is a script, when it does not start with the term "function". So simply omit the "format longG " or move it inside the function.

Notes:

• 1.79*10^-7 is an expensive power operation, while 1.79E-7 is a cheap constant.
• Using a pile of global variables whose names contain numbers is a really bad idea. It will be nmearly impossible to debug this. Prefer using indices and indices instead of hiding indices in the name of variables. The numbers variables lead to the really frightening code.

[Avani Dwivedi]

Thanks a lot Jan Simon. But please can you give a look to the file attached above. And could you please explain by - 'You cannot create functions inside a script or in the command window. An M-file is a script, when it does not start with the term "function"'. Can you give the example from the code? Thanking you in advance.