How to convert my script to a function

Hello everyone, thanks so much for all the help I am getting here to learn MATLAB. Please how can I convert my script into a function. I did an integration using Euler explicit method and I want to use it for the objective function of my optimization, so I want to convert the script into a function. Please I really need help with this.

%integration of gas-lift ODEs using euler explicit methods

% clear all

% clc

% close

%timesteps

tf = 3600;

h = 0.1;

t = 0:h:tf;

n = length(t);

%initialising

x10 = 1.0568;

x20 = 1.0606;

x30 = 0.6000;

x40 = 0.6700;

x50 = 6.0000;

x60 = 5.5047;

x70 = 0.1265;

x80 = 0.9863;

x(1,1) = x10;

x(2,1) = x20;

x(3,1) = x30;

x(4,1) = x40;

x(5,1) = x50;

x(6,1) = x60;

x(7,1) = x70;

x(8,1) = x80;

%parameters

L_w = 1500;

H_w = 1000;

D_w = 0.121;

L_bh = 500;

H_bh = 500;

D_bh = 0.121;

L_a = L_w;

H_a = H_w;

D_a = 0.189;

V_a = L_a*pi*((D_a/2)^2 - (D_w/2)^2);

A_w = pi*(D_w/2)^2;%[m2]

A_bh = pi*(D_bh/2)^2;%[m2]

L_r = 500;

H_r = 500;

D_r = 0.121;

A_r = pi*(D_r/2)^2;

rho_o = 8*1e2;

rho_ro = 8*1e2;

C_iv = 0.1e-3;

C_pc = 2e-3;

C_rh = 10e-3;

GOR1 = 0.1;

GOR2 = 0.12;

dp1 = 1e-6;

dp2 = 1e-6;

p_res1 = 150;

p_res2 = 155;

PI = 0.7;

T_a = 28+273;

T_w = 32+273;

T_r = 30+273;

p_s = 20;

% p_m = p_s;

mu_o = 0.001;

slip = 1;

R = 8.3145;

Mw = 0.020;

g = 9.8;

w_gl1 =1;

w_gl2 =1;

H_k = [0 0 0 0 0 0 20642.89 -374.645 0 0; 0 0 0 0 0 0 3674.553 -163.461 0 0];

w_k = [0.001 0.001 0.001 0.001 0.02 0.02 0.01 0.03 0.001 0.001]' ;

v_k = rand(2,1);

q = cov(w_k);

r = cov(v_k);

R_k = [0.01 0; 0 0.01];

lb = [10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-2;10*1e-2];

ub = [10e7;10e7;10e7;10e7;10e7;10e7;10e7;10e7;50e4;50e4];

QgMax = 10;

qGLMax = 40;

a = 1;

b = 0.5;

A = [];

b = [];

Aeq = [];

beq = [];

for i=1:n

%auxiliary algebraic equations

p_a1 = ((R*T_a)/(V_a*Mw) + (g*H_a/V_a))*x10*1e3;%pressure of gas in the annulus

p_a2 = ((R*T_a)/(V_a*Mw) + (g*H_a/V_a))*x20*1e3;

rho_a1 = (Mw/(R*T_a))*p_a1;%density of gas in the annulus

rho_a2 = (Mw/(R*T_a))*p_a2;

p_wh1 = ((R*T_w/Mw)*(x30*1e3/(L_w*A_w + L_bh*A_bh - x50*1e3/rho_o))) - (((x30*1e3+x50*1e3)/(L_w*A_w))*g*H_w/2);%pressure at the wellhead

p_wh2 = ((R*T_w/Mw)*(x40*1e3/(L_w*A_w + L_bh*A_bh - x60*1e3/rho_o))) - (((x40*1e3+x60*1e3)/(L_w*A_w))*g*H_w/2);

rho_m1 = (((x30*1e3+x50*1e3)*p_wh1*Mw*rho_o)/(x50*1e3*p_wh1*Mw + rho_o*R*T_w*x30*1e3));

rho_m2 = (((x40*1e3+x60*1e3)*p_wh2*Mw*rho_o)/(x60*1e3*p_wh2*Mw + rho_o*R*T_w*x40*1e3));

p_rh = (R*T_r/Mw)*(x70*1e3/(L_r*A_r)) - (((x70*1e3+x80*1e3)/(L_r*A_r))*g*H_r/2); %25 Pressure at the riser head 1[pa]

rho_r = ((x70*1e3+x80*1e3)*p_rh*Mw*rho_ro)/(x80*1e3*p_rh*Mw + rho_ro*R*T_r*x70*1e3); % density of mixture at the riser [pa]

w_pr = C_rh*sqrt(rho_r*(p_rh - p_s*1e5)); %28 total mass flow rate through the production choke1[kg/s]

p_m = p_rh + (g/(A_r*L_r)*(x80*1e3+x70*1e3)*H_r) + (128*mu_o*L_r*w_pr/(pi*D_r^4*((x70*1e3 + x80*1e3)*p_rh*Mw*rho_ro)/(x80*1e3*p_rh*Mw + rho_ro*R*T_r*x70*1e3)));%27%Pressure at the manifold[pa]

w_pc1 = C_pc*sqrt(rho_m1*(max(1e-3,(p_wh1 - p_m))));

w_pc2 = C_pc*sqrt(rho_m2*(max(1e-3,(p_wh2 - p_m))));

w_pg1 = (x30*1e3/(max(1e-3,(x30*1e3+x50*1e3))))*w_pc1;

w_pg2 = (x40*1e3/(max(1e-3,(x40*1e3+x60*1e3))))*w_pc2;

w_po1 = (x50*1e3/(max(1e-3,(x30*1e3+x50*1e3))))*w_pc1;

w_po2 = (x60*1e3/(max(1e-3,(x40*1e3+x60*1e3))))*w_pc2;

p_w1 = p_wh1 +((g/(A_w*L_w))*(max(1e-3,(x30*1e3+x50*1e3-rho_o*L_bh*A_bh)))*H_w)+ (128*mu_o*L_w*w_pc1/(pi*D_w^4*(x30*1e3+x50*1e3)*p_wh1*Mw*rho_o)/(x50*1e3*p_wh1*Mw + rho_o*R*T_w*x30*1e3)); %((m_gt1 + m_ot1)*p_wh1*Mw.*rho_o)./(m_ot.*1e3.*p_wh.*1e5.*Mw + rho_o.*R.*T_w.*m_gt.*1e3)));;

p_w2 = p_wh2 +((g/(A_w*L_w))*(max(1e-3,(x40*1e3+x60*1e3-rho_o*L_bh*A_bh)))*H_w)+ (128*mu_o*L_w*w_pc2/(pi*D_w^4*(x40*1e3+x60*1e3)*p_wh2*Mw*rho_o)/(x60*1e3*p_wh2*Mw + rho_o*R*T_w*x40*1e3));

w_iv1 = C_iv*sqrt(rho_a1*(max(1e-3,(p_a1 - p_w1))));

w_iv2 = C_iv*sqrt(rho_a2*(max(1e-3,(p_a2 - p_w2))));

p_bh1 = p_w1 + (rho_o*g*H_bh) + (128*mu_o*L_bh*w_po1/(pi*D_bh^4*rho_o));

p_bh2 = p_w2 + (rho_o*g*H_bh) + (128*mu_o*L_bh*w_po2/(pi*D_bh^4*rho_o));

w_ro1 = PI*1e-5*(p_res1*1e5 - p_bh1);

w_ro2 = PI*1e-5*(p_res2*1e5 - p_bh2);

w_rg1 = GOR1*w_ro1;

w_rg2 = GOR2*w_ro2;

w_to = (x80*1e3/(x70*1e3+x80*1e3))*w_pr; % 29 mass flow rate of oil from the reservoir2 [kg/s]

w_tg = (x70*1e3/(x70*1e3+x80*1e3))*w_pr; %30 mass flow rate of gas from the reservoir [kg/s]

x(1,i+1) = x(1,i) + h*(w_gl1 - w_iv1)*1e-3;

x(2,i+1) = x(2,i) + h*(w_gl2 - w_iv2)*1e-3;

x(3,i+1) = x(3,i) + h*(w_iv1 + w_rg1 - w_pg1)*1e-3;

x(4,i+1) = x(4,i) + h*(w_iv2 + w_rg2 - w_pg2)*1e-3;

x(5,i+1) = x(5,i) + h*(w_ro1 - w_po1)*1e-3;

x(6,i+1) = x(6,i) + h*(w_ro2 - w_po2)*1e-3;

x(7,i+1) = x(7,i) + h*((w_pg1+w_pg2)-w_tg)*1e-3;

x(8,i+1) = x(8,i) + h*((w_po1+w_po2)-w_to)*1e-3;

x_hat_aug = [x(1,i+1);x(2,i+1);x(3,i+1);x(4,i+1);x(5,i+1);x(6,i+1);x(7,i+1);x(8,i+1);GOR1;GOR2];

h_k = (H_k*x_hat_aug);

P_0 = diag(ones(size(x_hat_aug)));

y_meas = h_k + 1e-3*h_k*randn;

P_k =F_k*P_0*F_k'+Q;

K_k = P_k*H_k'*inv(H_k*P_k*H_k'+R_k);

x_bark = x_hat_aug + K_k*(y_meas - h_k);

l=eye(10);

P_kk = (l-K_k*H_k)*P_k;

x_est =x_bark;

X_est_all(i,:)=x_est;

y_est_all(i,:)=y_meas;

h_k_all(i,:)= h_k;

% % y1 = w_to(:,end);

% y2 = w_tg(:,end);

p = [x_est(9) x_est(10)];

all_p(i,:)= p;

x_aug = [x(1,i);x(2,i);x(3,i);x(4,i);x(5,i);x(6,i);x(7,i);x(8,i);GOR1;GOR2];

x10 = x(1,i+1);

x20 = x(2,i+1);

x30 = x(3,i+1);

x40 = x(4,i+1);

x50 = x(5,i+1);

x60 = x(6,i+1);

x70 = x(7,i+1);

x80 = x(8,i+1);

[Adam]

Just put

function myFunction( )

at the top, replacing 'myfunction' there with whatever you called the .m file.

If you need to pass in arguments then these would go between the ( and ), but that code is far too long for me to want to try to work out if it wants input arguments or not!

The documentation also gives help on using functions as they are a fundamental part of the language.