title: Cubic Spline Program
body: <p>I'm trying to write a cubic spline interpolation program. I have written the program but, the graph is not coming out correctly. The spline uses natural boundary conditions(second dervative at start/end node are 0). The code is in Matlab and is shown below,</p>
  
  <pre><code>clear all
  %Function to Interpolate
  k = 10;                    %Number of Support Nodes-1
  xs(1) = -1;
  for j = 1:k
      xs(j+1) = -1 +2*j/k;   %Support Nodes(Equidistant)
  end;
  fs = 1./(25.*xs.^2+1);     %Support Ordinates
  x = [-0.99:2/(2*k):0.99];  %Places to Evaluate Function
  fx = 1./(25.*x.^2+1);      %Function Evaluated at x
  
  %Cubic Spline Code(Coefficients to Calculate 2nd Derivatives)
  
  f(1) = 2*(xs(3)-xs(1));
  g(1) = xs(3)-xs(2);
  r(1) = (6/(xs(3)-xs(2)))*(fs(3)-fs(2)) + (6/(xs(2)-xs(1)))*(fs(1)-fs(2));
  e(1) = 0;
  
  for i = 2:k-2
      e(i) = xs(i+1)-xs(i);
      f(i) = 2*(xs(i+2)-xs(i));
      g(i) = xs(i+2)-xs(i+1);
      r(i) = (6/(xs(i+2)-xs(i+1)))*(fs(i+2)-fs(i+1)) + ...
             (6/(xs(i+1)-xs(i)))*(fs(i)-fs(i+1));
  end
  e(k-1) = xs(k)-xs(k-1);
  f(k-1) = 2*(xs(k+1)-xs(k-1));
  r(k-1) = (6/(xs(k+1)-xs(k)))*(fs(k+1)-fs(k)) + ...
           (6/(xs(k)-xs(k-1)))*(fs(k-1)-fs(k));
  
  %Tridiagonal System
  
  i = 1;
  A = zeros(k-1,k-1);
  while i &lt; size(A)+1;
      A(i,i) = f(i);
      if i &lt; size(A);
          A(i,i+1) = g(i);
          A(i+1,i) = e(i);
      end
      i = i+1;
  end
  
  for i = 2:k-1                             %Decomposition
      e(i) = e(i)/f(i-1);
      f(i) = f(i)-e(i)*g(i-1);
  end
  
  for i = 2:k-1                             %Forward Substitution 
      r(i) = r(i)-e(i)*r(i-1);
  end
  
  xn(k-1)= r(k-1)/f(k-1);
  for i = k-2:-1:1                          %Back Substitution
      xn(i) = (r(i)-g(i)*xn(i+1))/f(i);
  end
  
  %Interpolation
  
  if (max(xs) &lt;= max(x))
      error('Outside Range'); 
  end
  
  if (min(xs) &gt;= min(x))
      error('Outside Range'); 
  end
  
  
  P = zeros(size(length(x),length(x)));
  i = 1;
  for Counter = 1:length(x)
      for j = 1:k-1
          a(j) = x(Counter)- xs(j);
      end
      i = find(a == min(a(a&gt;=0)));
      if i == 1
          c1 = 0;
          c2 = xn(1)/6/(xs(2)-xs(1));
          c3 = fs(1)/(xs(2)-xs(1));
          c4 = fs(2)/(xs(2)-xs(1))-xn(1)*(xs(2)-xs(1))/6;
          t1 = c1*(xs(2)-x(Counter))^3;
          t2 = c2*(x(Counter)-xs(1))^3;
          t3 = c3*(xs(2)-x(Counter));
          t4 = c4*(x(Counter)-xs(1));
          P(Counter) = t1 +t2 +t3 +t4;
      else
          if i &lt; k-1
          c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
          c2 = xn(i+1)/6/(xs(i+1)-xs(i-1+1));
          c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
          c4 = fs(i+1)/(xs(i+1)-xs(i-1+1))-xn(i+1)*(xs(i+1)-xs(i-1+1))/6;
          t1 = c1*(xs(i+1)-x(Counter))^3;
          t2 = c2*(x(Counter)-xs(i-1+1))^3;
          t3 = c3*(xs(i+1)-x(Counter));
          t4 = c4*(x(Counter)-xs(i-1+1));
          P(Counter) = t1 +t2 +t3 +t4;
          else
          c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
          c2 = 0;
          c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
          c4 = fs(i+1)/(xs(i+1)-xs(i-1+1));
          t1 = c1*(xs(i+1)-x(Counter))^3;
          t2 = c2*(x(Counter)-xs(i-1+1))^3;
          t3 = c3*(xs(i+1)-x(Counter));
          t4 = c4*(x(Counter)-xs(i-1+1));
          P(Counter) = t1 +t2 +t3 +t4;
          end    
      end
  end
  
  P = P';
  P(length(x)) = NaN;
  
  plot(x,P,x,fx)
  </code></pre>
  
  <p>When I run the code, the interpolation function is not symmetric and, it doesn't converge correctly. Can anyone offer any suggestions about problems in my code? Thanks.</p>