How to repeat a script for variations of a parameter

Hello everyone, Im writting a script for the blade element momentum theory as follows:

if true
  clc;
clear all;
load('H:\tmp\Matlab\datos con tabla')
clear Cl Cd
  airdensity = 1.25;
  pitchangle = 0;
  B = 3;                    %number of blades
  pi = 3.1416;
  revolutions=17.5;
  omega = (17.5*2*pi)/60;
  R = 40;                   %radius of the rotor
  Vo=11.9
  Nelements = 23;
  bladetwist=bladetwist(:,1)
  r = [1.24 3.24 5.24 7.24 9.24 11.24 13.24 15.24 17.24 19.24 21.24 23.24 25.24 27.24 29.24 31.24 33.24 35.24 37.24 38.24 39.24 39.64 40.04];  
  c = [0.1 0.5 1 2.81 2.98 3.14 3.17 2.99 2.79 2.58 2.38 2.21 2.06 1.92 1.8 1.68 1.55 1.41 1.18 0.98 0.62 0.48 0.07];
  dy = [1.24 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 0.4 0.4];
  Vtip = omega*R;
  tsr = Vtip/Vo;             %tip speerd ratio
  n = length(r);
Ptot = 0;
Mtot = 0;
Ttot = 0;
%Cn = zeros(1,n);
%Ct = zeros(1,n);
axif = zeros(1,n);
tanif = zeros(1,n);
%dPn = zeros(1,n);
%dPt = zeros(1,n);
Mtot = zeros(1,n);
Ttot = zeros(1,n);
Ptot = zeros(1,n);
%dM = zeros(1,n);
%dT = zeros(1,n);
%dFl = zeros(1,n);
%dFd = zeros(1,n);
Fltot = 0;
Fdtot = 0;
isconverged=0;
axif=zeros(1,n);
tanif=zeros(1,n);
newaxif=zeros(1,n);
newtanif=zeros(1,n);
tolerance=0.000001;
for i = 1:n
 while ~ isconverged
  inflowangle(i)=atan((Vo(s)*(1-newaxif(i)))/(omega*r(i)*(1+newtanif(i))));
  inflowangle_degrees(i)=inflowangle(i)*(360/(2*pi));
  AOA_degrees(i)=inflowangle_degrees(i)-pitchangle-bladetwist(i);
  AOA(i)=inflowangle(i)-pitchangle;

    x=polarCl(:,1);
    y=polarCl(:,i+1);
    xo=AOA_degrees(i);
    Cl(i)=interp1(x,y,xo);

    x1=polarCd(:,1);
    y1=polarCd(:,i+1);
    xo1=AOA_degrees(i);
    Cd(i)=interp1(x1,y1,xo1);

    Cn(i)=Cl(i)*cos(inflowangle(i))+Cd(i)*sin(inflowangle(i));
    Ct(i)=Cl(i)*sin(inflowangle(i))-Cd(i)*cos(inflowangle(i));

    sigma(i)=(c(i)*B)/(2*pi*r(i));
    axif(i)=newaxif(i);
    tanif(i)=newtanif(i);
    newaxif(i)=1/(((4*sin(inflowangle(i))^2)/(sigma(i)*Cn(i)))+1);
    newtanif(i)=1/(((4*sin(inflowangle(i))*cos(inflowangle(i)))/(sigma(i)*Ct(i)))-1);
    if (abs(newaxif(i)-axif(i))<tolerance) && (abs(newtanif(i)-tanif(i))<tolerance)
        isconverged=1;
    end

   end
   isconverged=0;

    Vabs(i)=Vo(s)/sin(inflowangle(i));

    dFl(i)=0.5*airdensity*Vabs(i)^2*c(i)*Cl(i)*dy(i)
    Fltot=Fltot+dFl(i);

    dFd(i)=0.5*airdensity*Vabs(i)^2*c(i)*Cd(i)*dy(i)
    Fdtot=Fdtot+dFd(i);

    dPn(i)=dFl(i)*cos(inflowangle(i))+dFd(i)*sin(inflowangle(i));
    dPt(i)=dFl(i)*sin(inflowangle(i))-dFd(i)*cos(inflowangle(i));

    f(i)=(B*(R-r(i)))/(2*r(i)*sin(inflowangle(i)));
    F(i)=(2*cos(exp(-f(i)))^-1)/pi;

    dT(i)=B*dPn(i)*F(i)
    Ttot = Ttot + dT(i);

    dM(i)=r(i)*B*dPt(i)*F(i)
    Mtot = Mtot + dM(i);

    dP(i)=dM(i)*omega
    Ptot = Ptot + dP(i);
end

figure (1) plot(r,dT) xlabel('r') ylabel('dT') figure (2) plot(r,dM) xlabel('r') ylabel('dM') figure (3) plot(r,dP) xlabel('r') ylabel('dP') end

now i would like to write another script to repeat this last one many times for diferent values of Vo, and i wrote this but its not working:

if true
clc
close all
clear all
load('H:\tmp\Matlab\datos con tabla')
clear Cl Cd
Vo=[3:0.5:18];
z=length(Vo);
Ptot=zeros(1,z);
Mtot=zeros(1,z);
Ttot=zeros(1,z);
for s = 1:z

   [Ptot(s), Mtot(s), Ttot(s)] = pruevaliteratura_convergence(Vo(s))
end
  end

It tells me the next error: Attempt to execute SCRIPT pruevaliteratura_convergence as a function: H:\tmp\Matlab\pruevaliteratura_convergence.m

Error in curva_potencia (line 13) [Ptot(s), Mtot(s), Ttot(s)] = pruevaliteratura_convergence(Vo(s))

I dont know why its not working, I just want get the different outputs for diferent Vo and then plot them as a function of Vo.

Thank you beforehand.

[Dennis]

I think the cleanest approach would be to write a function instead of using a script.

[Enrique Montero Mora]

but if i want to use the for loop in the first script, dont i need to write it as an script?

[Stephen23]

I agree with Dennis: the best solution would be to abandon scripts (which is the best thing to do with scripts anyway) and convert the code to a function. Then your task is trivial: call the function in a loop, collect the outputs, plot.

[Enrique Montero Mora]

and any idea of how can I do it? Im new in matlab and i have no clue of using functions...

[Torsten]

function main
  Vo=[3:0.5:18];
  z=length(Vo);
  Ptot=zeros(1,z);
  Mtot=zeros(1,z);
  Ttot=zeros(1,z);
  for s = 1:z
    [Ptot(s), Mtot(s), Ttot(s)] = pruevaliteratura_convergence(Vo(s))
  end
end

function [Ptot, Mtot, Ttot] = pruevaliteratura_convergence(Vo)
  1. Insert your code from above 
  2. Remove the line "Vo=11.9"
  3. Replace
    Vabs(i)=Vo(s)/sin(inflowangle(i));
  by
    Vabs(i)=Vo/sin(inflowangle(i)); 
  4. If possible, don't repeatedly load 'H:\tmp\Matlab\datos con tabla')
end

[Stephen23]

https://www.mathworks.com/help/matlab/matlab_prog/create-functions-in-files.html

[dpb]

Others have commented on and shown how to write functions vis a vis scripts; I'll just note that

pi = 3.1416;

overloads the Matlab builtin in pi function that returns the value of pi to full double precision with a crude approximation of not even single precision accuracy. Delete the above line and use pi where want the value.

[Enrique Montero Mora]

Thank you very much guys! great help! :) @torsten I tried to do what you told me, but up the fifth iteration it wont stop debugging the interpol1 function and gets stucked there. I did this:

function main
  Vo=[3:0.5:14];
  z=length(Vo);
  Ptot=zeros(1,z);
  Mtot=zeros(1,z);
  Ttot=zeros(1,z);
  for s = 1:z
    [Ptot(s), Mtot(s), Ttot(s)] = pruevaliteratura_convergence(Vo(s))
  end
end

function [Ptot, Mtot, Ttot] = pruevaliteratura_convergence(Vo)
  airdensity = 1.25;
  pitchangle = 0;
  B = 3;                    %number of blades
  revolutions=17.5;
  omega = (17.5*2*pi)/60;
  R = 40;                   %radius of the rotor
  Nelements = 23;
  bladetwist=bladetwist(:,1)
  r = [1.24 3.24 5.24 7.24 9.24 11.24 13.24 15.24 17.24 19.24 21.24 23.24 25.24 27.24 29.24 31.24 33.24 35.24 37.24 38.24 39.24 39.64 40.04];  
  c = [0.1 0.5 1 2.81 2.98 3.14 3.17 2.99 2.79 2.58 2.38 2.21 2.06 1.92 1.8 1.68 1.55 1.41 1.18 0.98 0.62 0.48 0.07];
  dy = [1.24 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 0.4 0.4];
  Vtip = omega*R;
  %tip speerd ratio
  n = length(r);
  Ptot = 0;
  Mtot = 0;
  Ttot = 0;
  %Cn = zeros(1,n);
  %Ct = zeros(1,n);
  %dPn = zeros(1,n);
  %dPt = zeros(1,n);
  %Mtot = zeros(1,n);
  %Ttot = zeros(1,n);
  %Ptot = zeros(1,n);
  %dM = zeros(1,n);
  %dT = zeros(1,n);
  %dFl = zeros(1,n);
  %dFd = zeros(1,n);
  Fltot = 0;
  Fdtot = 0;
  isconverged=0;
  axif=zeros(1,n);
  tanif=zeros(1,n);
  newaxif=zeros(1,n);
  newtanif=zeros(1,n);
  tolerance=0.000001;
  for i = 1:n
    while ~ isconverged
      inflowangle(i)=atan((Vo*(1-newaxif(i)))/(omega*r(i)*(1+newtanif(i))));
      inflowangle_degrees(i)=inflowangle(i)*(360/(2*pi));
      AOA_degrees(i)=inflowangle_degrees(i)-pitchangle-bladetwist(i);
      AOA(i)=inflowangle(i)-pitchangle;
      x=polarCl(:,1);
      y=polarCl(:,i+1);
      xo=AOA_degrees(i);
      Cl(i)=interp1(x,y,xo);
      x1=polarCd(:,1);
      y1=polarCd(:,i+1);
      xo1=AOA_degrees(i);
      Cd(i)=interp1(x1,y1,xo1);

    Cn(i)=Cl(i)*cos(inflowangle(i))+Cd(i)*sin(inflowangle(i));
    Ct(i)=Cl(i)*sin(inflowangle(i))-Cd(i)*cos(inflowangle(i));

    sigma(i)=(c(i)*B)/(2*pi*r(i));
    axif(i)=newaxif(i);
    tanif(i)=newtanif(i);
    newaxif(i)=1/(((4*sin(inflowangle(i))^2)/(sigma(i)*Cn(i)))+1);
    newtanif(i)=1/(((4*sin(inflowangle(i))*cos(inflowangle(i)))/(sigma(i)*Ct(i)))-1);
    if (abs(newaxif(i)-axif(i))<tolerance) && (abs(newtanif(i)-tanif(i))<tolerance)
        isconverged=1;
    end

   end
   isconverged=0;

    Vabs(i)=Vo/sin(inflowangle(i));

    dFl(i)=0.5*airdensity*Vabs(i)^2*c(i)*Cl(i)*dy(i)
    Fltot=Fltot+dFl(i);

    dFd(i)=0.5*airdensity*Vabs(i)^2*c(i)*Cd(i)*dy(i)
    Fdtot=Fdtot+dFd(i);

    dPn(i)=dFl(i)*cos(inflowangle(i))+dFd(i)*sin(inflowangle(i));
    dPt(i)=dFl(i)*sin(inflowangle(i))-dFd(i)*cos(inflowangle(i));

    f(i)=(B*(R-r(i)))/(2*r(i)*sin(inflowangle(i)));
    F(i)=(2*cos(exp(-f(i)))^-1)/pi;

    dT(i)=B*dPn(i)*F(i)
    Ttot = Ttot + dT(i);

    dM(i)=r(i)*B*dPt(i)*F(i)
    Mtot = Mtot + dM(i);

    dP(i)=dM(i)*omega
    Ptot = Ptot + dP(i);
end

figure (1) plot(r,dT) xlabel('r') ylabel('dT') figure (2) plot(r,dM) xlabel('r') ylabel('dM') figure (3) plot(r,dP) xlabel('r') ylabel('dP') end end