function  [SP1,SP2,SP3] = fsplineiter(spa,spb,spc,x,y,knots,M0,kin,s,lambda1,AAe,bbe,AAne,bbne)

%Start
splines{1} = spa;
splines{2} = spb;
splines{3} = spc;
n = length(knots);
m = size(x,1);
ns = length(s);

%--------------------------------------------------------
%Destination vector
%--------------------------------------------------------
Y = [ y(:,1); y(:,2); y(:,3); ones(ns,1)*M0*lambda1 ];
X = zeros(m+m+m+ns,2*n+2*n+2*n);

%--------------------------------------------------------
%First part of X matrix: A,B,C points
%--------------------------------------------------------
for k=1:3,
  j = 1;
  while (j<=m),
    %search index of interval
    i = 1;
    while (i<n) & (x(j,k)>knots(i+1)),
      i = i+1;
    end
    if (i>=n) | (x(j,k)<knots(i)),
      error('Error: Out of intervalls!');
    end
    %calculate linear parameters
    k1 = knots(i); k2 = knots(i+1); xx = x(j,k); h = k2-k1;
    a = (k2-xx)^2*(xx-k1)/h^2; b = -(xx-k1)^2*(k2-xx)/h^2;
    c = (k2-xx)^2*(2*(xx-k1)+h)/h^3; d = (xx-k1)^2*(2*(k2-xx)+h)/h^3;
    %pack them into X matrix
    X((k-1)*m+j,(k-1)*2*n+i*2-1) = c;
    X((k-1)*m+j,(k-1)*2*n+i*2+0) = a;
    X((k-1)*m+j,(k-1)*2*n+i*2+1) = d;
    X((k-1)*m+j,(k-1)*2*n+i*2+2) = b;
    % next point in x()
    j = j+1;
  end
end

%--------------------------------------------------------
% Second part of X matrix: ca+..+k1*Ica+... = M0
%--------------------------------------------------------
j = 1;
while (j<=ns),
  for k=1:3,
    %search index of interval
    i = 1;
    while (i<n) & (s(j)>knots(i+1)),
      i = i+1;
    end
    if (i>=n) | (s(j)<knots(i)),
      error('Error: Out of intervalls!');
    end
    %calculate linear parameters
    k1 = knots(i); k2 = knots(i+1); xx = s(j); h = k2-k1;
    a = (k2-xx)^2*(xx-k1)/h^2; b = -(xx-k1)^2*(k2-xx)/h^2;
    c = (k2-xx)^2*(2*(xx-k1)+h)/h^3; d = (xx-k1)^2*(2*(k2-xx)+h)/h^3;
    [ee,aa,bb,cc,dd] = intspline2(splines{k},xx);
    %add to X
    rw = 3*m+j;  cl = (k-1)*2*n+i*2-1;
    X(rw,cl+0) = X(rw,cl+0) + lambda1*(c+kin(k)*cc);
    X(rw,cl+1) = X(rw,cl+1) + lambda1*(a+kin(k)*aa);
    X(rw,cl+2) = X(rw,cl+2) + lambda1*(d+kin(k)*dd);
    X(rw,cl+3) = X(rw,cl+3) + lambda1*(b+kin(k)*bb);
    Y(rw) = Y(rw) - lambda1*kin(k)*ee;
  end
  % next point in s()
  j = j+1;
end

%--------------------------------------------------------
% Solve least-squares problem
%--------------------------------------------------------

T = lsqlin(X,Y,AAne,bbne,AAe,bbe);

%Spline - A
SP1.x = knots;
SP1.y = T([1:2:2*n-1]);
SP1.dy = T([2:2:2*n]);
%Spline - B
SP2.x = knots;
SP2.y = T([2*n+1:2:4*n-1]);
SP2.dy = T([2*n+2:2:4*n]);
%Spline - C
SP3.x = knots;
SP3.y = T([4*n+1:2:6*n-1]);
SP3.dy = T([4*n+2:2:6*n]);