Help with my non-linear numeric code lift line method

Hello I have been working in this code based on the steps of the theory according to the aerodynamics book anderson

I have separately created the codes for the functions.

However, I don't know how to perform step 3, where I have to evaluate the induced angle for each station whit Numerical Simpson Rule

Can someone tell me how can I do it ,please?

%%Code Nonlinear Lifting Line:

clear;

clc;

format shorte;

b =10 ; %envergadura

cr =1.0 ;

ct =1.0 ;

twist =0.0 ;

N =10 ;%Numero de divisiones

rho =1.225 ; %densidad del aire

mu =1.7894e-5 ;

V =25.0 ;

%alpha angulo de ataque del ala

alpha =2.0 ;

%D factor de amortiguamiento

D=0.05 ;

[S,AR,lambda,y,c,alphag] = WingGeometry( b , cr , ct , twist , N);

alphainc=alphag+alpha ;

Rec = rho*V*c/mu ;

Index=1:N+1;

display(' Estacion y c alphag alphainc Rec');

display(' ------------------------------------------------------------------------------------------------------');

[Index' y' c' alphag' alphainc' Rec']

Gmax = 3.0 ;

GammaInput = sqrt(Gmax^2.0 *(1.0 - y.^2.0 ./ (b/2.0)^2.0));

GammaOld = GammaInput ;

dgdy = dydxFD(y ,GammaOld) ;

alphaind = EvaluateAlphaInd( y ,dgdy , V ) ;

alphaeff = alphainc - alphaind ;

PolarFileRead;

clr = interpl(alphaf,clf,alphaeff,'spline')

GammaNew = 0.5 * V .* c .* clr

GammaNew(1) = 0.0 ;

GammaNew(end) = 0.0 ;

GammaInput=GammaOld+D*(GammaNew-GammaOld);

%%%Code DydxFD :

function dydx = dydxFD(x,fx)

n = length(x);

dydx(1:n)-0.0;

dx = x(2) - x(1);

dydx(1) = (-fx(3) + 4.0*fx(2) -3.0*fx(1))/ (2.0 * dx );

for i=2:n-1

dydx(i) = ( fx(i+1) - fx(i-1) ) / (2.0 * dx );

end

dydx(n) - ( fx(n-2)- 4.0*fx(n-1)+ 3.0*fx(n))/ (2.0 * dx );

end

%%%Code Simpson_3:

function Integral = Simpson1_3( x , fx )

dx = x(2) - x(1) ;

n = lenght(x) ;

if mod(n-1,2) --0

Integral = 0 ;

for i=2:2:n-1

Integral = Integral + fx(l-1) + 4.0 * fx(i) + fx(l+1) ;

end

Integral= Integral * dx / 3.0 ;

else

display('Error: numero de puntos es par ! ') ;

Integral = 0.0 ;

end

end

%%%WingGeometry

function [S,AR,lambda,y,c,alphag] = WingGeometry( b , cr , ct , twist, N)

S=b * ( cr + ct ) /2.0 ;

AR = b ^2.0 / S ;

lambda = ct / cr ;

y=-b/2.0 : b/N : b/2 ;

c(1:N+1) = ct ;

c(N/2+1) = cr ;

c(2:N)= cr - (cr - ct)*abs(y(2:N))./(b/2.0) ;

alphag=twist*abs(y)./(b/2.0);

end

%%%Code PolarFile:

PolarFile= 'naca0012';

InputFile=fopen('naca0012');

for l=1:12;

DummyLine-fgetl(InputFile);

end

Read = 1 ;

Index=0;

while Read

Line-fgetl(InputFile);

if Line > 0

Index=Index+1;

Vars=sscanf(Line,'%If');

alphaf(Index)-Vars(1);

clf(Index)-Vars(2);

cdf(Index)-Vars(3);

cmf(Index)=Vars(5);

else

Read=0;

end

end

fclose(InputFile)

[Om prakash Meena]

i am also facing problem in it

[Om prakash Meena]

if someone have solutation for this problem, please reply on ""[e-mail]""