How to interpolate around an aerofoil.

Hi I have a piece of code which must interpolate/loft aerofoil sections together. Each of these sections may have it's own aerofoil so the points are scattred. I try interpolating and then scaling to get to the correct size but cannot get the leading edge of the aerofoil to interpolate properly. I'm using Triscattredinterp and seperating the top and bottom surfaces. I Have tried changing the limits for interpolation, and changing the interpolation method.

Below is my code as well as what I run it with.

% aerofoilcoord = getBladeGeometry(R,chords,aerofoils)
% Interpolate between the aerofoils then apply the chord scaling to each
% section.
% aerofoilcoord = the coordinates defining the blade geometry
% R = Length of the blade section of the prop, ie, R-RHub
% chords = The chords
% aerofoils = cell arraywith the xyz that defines each aerofoil
% npoints = Half the number of points to use on a section

function  aerofoilcoord = getBladeGeometry(R,chords,aerofoils,npoints)
% Thrust axis - y aerofoils
% Chord axis - z aerofoils
% span axis - x aerofoils
% Essential problem is there is a top and bottom value for each rpos cpos
% so split into top and bottom surfaces
if(nargin<4)
    npoints = length(aerofoils{1});
end
UpperSurface = [];
LowerSurface = [];
for a = 1:length(aerofoils)
    % find where chord changes direction
    V = [0 ;diff(aerofoils{a}(2:end,1))>0];
    if (V(1)==0)
        mpoint = find(V,1,'first');  % First index where V is +
    else
        mpoint = find(V,0,'first');  % First index where V is -
    end    
    % add aerofoil sections together
    UpperSurface = [UpperSurface; aerofoils{a}(1:mpoint,:)];
    LowerSurface = [LowerSurface; aerofoils{a}(mpoint:end,:)];
end

% r then chord, then height
U = TriScatteredInterp(UpperSurface(:,3),UpperSurface(:,1),UpperSurface(:,2),'natural');
L = TriScatteredInterp(LowerSurface(:,3),LowerSurface(:,1),LowerSurface(:,2),'natural');

r = linspace(R,0,length(chords)); % Want the geometry in tip to root order
%aerofoilcoord = zeros(length(chords)*2*npoints,3);
cpos = linspace(0,1,npoints);
yu = zeros(npoints,1);
yl = zeros(npoints,1);

% For Testing purposes I interpolate a single section here.
plot(cpos,U(ones(1,length(cpos)),cpos),'r')
hold on
plot(cpos,L(ones(1,length(cpos)),cpos),'b')

pos = 1;
for i = 1:length(chords)
    for j = 1:npoints
        yu(j) = U(r(i),cpos(j)); %Upper part
        yl(j) = L(r(i),cpos(j)); %Lower part
        if( isnan(yu(j))) % Doesn't like interpolating near to cpos =1
            yu(j) = yu(j-1);
        end
        if( isnan(yl(j)))
            yl(j) = yl(j-1);
        end
    end
    y = [flipud(yu) ; yl]; % maybe remove duplicate entries
    aerofoilcoord(pos:pos+npoints*2-1,1) = r(i);
    aerofoilcoord(pos:pos+npoints*2-1,2) = y*chords(i);
    aerofoilcoord(pos:pos+npoints*2-1,3) = [fliplr(cpos)  cpos]'*chords(i);

      pos = pos+npoints;
  end

end

chords = [0 0.0356700514213020 0.0508708806851791 0.0628139407496790 0.0731035002578730 0.0823489045764129 0.0908541115870836 0.0987926357852838 0.106271294763325 0.113358566428082 0.120098971323548 0.126521096371485 0.132642418712841 0.138472393969863 0.144014547873467 0.149267969763741 0.154228434406592 0.158889286040854 0.163242166070747 0.167277634580859 0.170986598432440 0.174393505481916 0.177523023333961 0.180389119458216 0.183006673166425 0.185391303971780 0.187559199392091 0.189526949743190 0.191311395565056 0.192929491702286 0.194398190709030 0.195734347153519 0.196954643544188 0.198075537970250 0.199113233124374 0.200083666132930 0.201002518539856 0.201885245855795 0.202747126280689 0.203603328526789 0.204457272681934 0.205210248756096 0.205835359440329 0.206330937493283 0.206695164847081 0.206926081784407 0.207021596236942 0.206979493369704 0.206797445614583 0.206473023317337 0.206003706164941 0.205386895564501 0.204619928150381 0.203700090602776 0.202624635968195 0.201390801680077 0.199995829485585 0.198436987492215 0.196711594554723 0.194817047228499 0.192750902589509 0.190511018145157 0.188095225915324 0.185501548254530 0.182728264911183 0.179773960101909 0.176637573639828 0.173318456246522 0.169816429129717 0.166131847847887 0.162265670407236 0.158219529444003 0.153995808234054 0.149597720140621 0.145029390958654 0.140295943439813 0.135403583085473 0.130359684077298 0.125172873977763 0.119853115580163 0.114385786798859 0.108738775297409 0.102940634436196 0.0970228446862634 0.0910178582794671 0.0849588746886982 0.0788796160300551 0.0728141027403287 0.0667964295159430 0.0608605412668931 0.0550400087584836 0.0493678036927000 0.0438760732180185 0.0385959142395852 0.0335571484094775 0.0287880992785582 0.0243153737485825 0.0201636506306400 0.0163554797444330 0.0129110955315890 ]

aerofoils = {Naca009,sd7037,Naca009} % Can be any aerofoils, I'll send the files if you want them or grab them from http://www.ae.illinois.edu/m-selig/ads/coord_database.html aerofoilcoords = getBladeGeometry(1.2,chords,aerofoils);

Hope I've been clear, please help any suggestions would be usefull. Thanks.