>    restart;
SecDeri := proc(x,y,h,n)
  local i,j,A,b;
  with(linalg):
  A:=array(1..n,1..n);
  b:=array(1..n);
  # Reset A to 0
  for i from 1 to n do
    for j from 1 to n do
      A[i,j]:=0;
    od;
  od;
  A[1,1]:=1; b[1]:=0; A[n,n]:=1; b[n]:=0;
  for i from 2 to n-1 do
    A[i,i-1]:=h[i-1]/6;
    A[i,i]:=(h[i-1]+h[i])/3;
    A[i,i+1]:=h[i]/6;
    b[i]:=(y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1];
  od;
  RETURN(linsolve(A,b));
end:

CubicSpline := proc()
  local fp,filename, N1,X1,Y1, N2,X2,Y2, F1,F2,P1,P2,
        i,j, H1,M1, H2,M2, a,b,c,d, Udp,Ldp,U,L;
  with(plots):
  printf("***** Author: Mengnien Wu *************************************************\n");
  printf("* Draw 2 cubic splines for N1 and N2 data points respectively.            *\n");
  printf("* Data points are supposed obtained from a graph that can be decomposed   *\n");
  printf("* into 2 parts, upper and lower, graphs of some continuous functions, and *\n");
  printf("* N1 points from the upper part and N2 from the lower one are picked.     *\n");
  printf("* A file which contains                                                   *\n");
  printf("*              N1, X1[1],Y1[1] ... X1[N1],Y1[N1],                         *\n");
  printf("*              N2, X2[1],Y2[1] ... X2[N2],Y2[N2]                          *\n");
  printf("* will be read.                                                           *\n");
  printf("***************************************************************************\n");
  printf("\nEnter your filename:\n");
  filename := scanf("%a")[1];
  fp:=fopen(filename,READ);
    N1:=fscanf(fp,"%d")[1];
    X1:=array(1..N1);
    Y1:=array(1..N1);
    for i from 1 to N1 do
      X1[i]:=fscanf(fp,"%f")[1];
      Y1[i]:=fscanf(fp,"%f")[1];
    od;
    N2:=fscanf(fp,"%d")[1];
    X2:=array(1..N2);
    Y2:=array(1..N2);
    for i from 1 to N2 do
      X2[i]:=fscanf(fp,"%f")[1];
      Y2[i]:=fscanf(fp,"%f")[1];
    od;
  fclose(fp);

  F1:=array(1..N1-1); #Function
  P1:=array(1..N1-1); #Picture
  H1:=array(1..N1-1);
  for i from 1 to N1-1 do
    H1[i]:=X1[i+1]-X1[i];
  od;
  M1:=SecDeri(X1,Y1,H1,N1);

  for j from 1 to N1-1 do
    a :=Y1[j];
    b :=(Y1[j+1]-Y1[j])/H1[j]-(2*M1[j]+M1[j+1])*H1[j]/6;
    c :=M1[j]/2;
    d :=(M1[j+1]-M1[j])/(6*H1[j]);
    F1[j]:=a+(x-X1[j])*( b+(x-X1[j])*( c+(x-X1[j])*d ) ); #Horner's
    P1[j]:=plot(F1[j],x=X1[j]..X1[j+1],color=blue,scaling=constrained):
  od;

  F2:=array(1..N2-1); #Function
  P2:=array(1..N2-1); #Picture
  H2:=array(1..N2-1);
  for i from 1 to N2-1 do
    H2[i]:=X2[i+1]-X2[i];
  od;
  M2:=SecDeri(X2,Y2,H2,N2);
  for j from 1 to N2-1 do
    a := Y2[j];
    b := (Y2[j+1]-Y2[j])/H2[j] - (2*M2[j]+M2[j+1])*H2[j]/6;
    c := M2[j]/2;
    d := (M2[j+1]-M2[j])/(6*H2[j]);
    F2[j]:= a+(x-X2[j])*( b+(x-X2[j])*( c+(x-X2[j])*d ) ); #Horner's
    P2[j]:= plot(F2[j],x=X2[j]..X2[j+1],color=blue,scaling=constrained):
  od;

  Udp:=plot([ [X1[n],Y1[n]] $n=1..N1 ],style=point,symbol=circle):
  Ldp:=plot([ [X2[n],Y2[n]] $n=1..N2 ],style=point,symbol=circle):
  display(Udp,Ldp,[P1[n] $n=1..N1-1],[P2[n] $n=1..N2-1]);
end:

>    CubicSpline();

Warning, the name changecoords has been redefined

***** Author: Mengnien Wu *************************************************

* Draw 2 cubic splines for N1 and N2 data points respectively.            *

* Data points are supposed obtained from a graph that can be decomposed   *

* into 2 parts, upper and lower, graphs of some continuous functions, and *

* N1 points from the upper part and N2 from the lower one are picked.     *

* A file which contains                                                   *

*              N1, X1[1],Y1[1] ... X1[N1],Y1[N1],                         *

*              N2, X2[1],Y2[1] ... X2[N2],Y2[N2]                          *

* will be read.                                                           *

***************************************************************************

Enter your filename:

>    "DataPoints.txt"

Warning, the protected names norm and trace have been redefined and unprotected

[Maple Plot]

>