1. c = 1;
  2. t = .12*c;
  3. delE = c/n;
  4. n = 2;
  5. Vinf = 1;
  6.  
  7. (*Place points on airfoil*)
  8. Xfoil = Table[(i - .99)*c/n, {i, n}]
  9. Yfoil = Table[(t/
  10.       0.2)*(0.2969*Sqrt[Xfoil[[i]]/c] - (0.126*Xfoil[[i]])/c - 
  11.       0.3516*(Xfoil[[i]]/c)^2 + 0.2843*(Xfoil[[i]]/c)^3 - 
  12.       0.1015*(Xfoil[[i]]/c)^4), {i, n}];
  13.  
  14. Ej = Table[(i - .5)*delE, {i, n}];
  15. amatrix = 
  16.   Table[Table[(Yfoil[[i]]*delE)/((Xfoil[[i]] - Ej[[j]])^2 + 
  17.        Yfoil[[i]]^2), {i, n}], {j, n}];
  18. xmatrix = LinearSolve[amatrix, (Yfoil*Vinf)];
  19. psi = Vinf*y - 
  20.    Sum[((y*delE)*xmatrix[[j]])/((x - Ej[[j]])^2 + y^2), {j, 1, n}]; 
  21. vi = D[psi, y]  ;
  22. vj = D[-psi, x]; 
  23. vxi = Table[[vi /. {x -> Xfoil[i], y -> Yfoil[i]}], {i, n}];
  24. vyi = Table[[vj /. {x -> Xfoil[i], y -> Yfoil[i]}], {i, n}];