c = 1;
t = .12*c;
delE = c/n;
n = 2;
Vinf = 1;
(*Place points on airfoil*)
Xfoil = Table[(i - .99)*c/n, {i, n}]
Yfoil = Table[(t/
0.2)*
(0.2969*
Sqrt[Xfoil
[[i
]]/c
] -
(0.126*Xfoil
[[i
]])/c -
0.3516*(Xfoil[[i]]/c)^2 + 0.2843*(Xfoil[[i]]/c)^3 -
0.1015*(Xfoil[[i]]/c)^4), {i, n}];
Ej = Table[(i - .5)*delE, {i, n}];
amatrix =
Table[Table[(Yfoil[[i]]*delE)/((Xfoil[[i]] - Ej[[j]])^2 +
Yfoil[[i]]^2), {i, n}], {j, n}];
xmatrix = LinearSolve[amatrix, (Yfoil*Vinf)];
psi = Vinf*y -
Sum[((y*delE)*xmatrix[[j]])/((x - Ej[[j]])^2 + y^2), {j, 1, n}];
vi = D[psi, y] ;
vj = D[-psi, x];
vxi = Table[[vi /. {x -> Xfoil[i], y -> Yfoil[i]}], {i, n}];
vyi = Table[[vj /. {x -> Xfoil[i], y -> Yfoil[i]}], {i, n}];