t=[0,1,2,3];
xx=linspace(-5,5);
plot(xx,BSplineEval(1,1,t,xx))
plot(xx,BSplineEval(1,1,t,xx),'o')
BSplineEval(1,1,t,1)

ans =

     0

plot(xx,BSplineEval(2,1,t,xx),'o')
plot(xx,BSplineEval(3,1,t,xx),'o')
plot(xx,BSplineEval(4,1,t,xx),'o')
??? Attempted to access t(5); index out of bounds because
numel(t)=4.

Error in ==> <a href="matlab: opentoline('/home/accounts/personale/clrmrc90/aa1213/advanced_numerical_analysis/BSplineEval.m',4,0)">BSplineEval at 4</a>
  yy(t(j) <= xx & xx < t(j+1)) = 1;


t

t =

     0     1     2     3

plot(xx,BSplineEval(3,1,t,xx),'o')
BSplineEval(3,1,t,3)

ans =

     1

BSplineEval(,1,t,3)
??? BSplineEval(,1,t,3)
                |
Error: Expression or statement is incorrect--possibly
unbalanced (, {, or [.

BSplineEval(1,1,t,3)

ans =

     0

BSplineEval(2,1,t,3)

ans =

     0

BSplineEval(3,1,t,3)

ans =

     0

BSplineEval(4,1,t,3)
??? Attempted to access t(5); index out of bounds because
numel(t)=4.

Error in ==> <a href="matlab: opentoline('/home/accounts/personale/clrmrc90/aa1213/advanced_numerical_analysis/BSplineEval.m',4,0)">BSplineEval at 4</a>
  yy(t(j) <= xx & xx < t(j+1)) = 1;

d = [-1,-0.5,0,0.6,1.4,0.7,-0.1,-0.8;...
0,0.3,0.5,0.4,0.1,-0.4,-0.5,-0.3];
plot(d(1,:),d(2,:),'*')
vertices=[1,2,3,4;5,6,7,8]

vertices =

     1     2     3     4
     5     6     7     8

repmat(vertices,2,1)

ans =

     1     2     3     4
     5     6     7     8
     1     2     3     4
     5     6     7     8

reshape(repmat(vertices,2,1),2,8)

ans =

     1     1     2     2     3     3     4     4
     5     5     6     6     7     7     8     8

curves2d
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves2d
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves2d
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves2d
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves2d
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves2d
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves2d
xx=linspace(-5,5);
plot(xx,N(1,xx-1))
plot(xx,N(1,xx-2))
plot(xx,N(1,xx-3))
plot(xx,N(1,xx-4))
xx=linspace(-10,10);
plot(xx,N(1,xx-5))
plot(xx,N(1,xx-6))
plot(xx,N(1,xx+3))
plot(xx,N(2,xx-1))
xx=linspace(-10,10,1000);
plot(xx,N(2,xx-1))
plot(xx,N(2,xx-2))
plot(xx,N(2,xx-3))
plot(xx,N(3,xx-1))
plot(xx,N(3,xx-3))
plot(xx,N(4,xx-1))
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves3d
d

d =

  Columns 1 through 5

         0    0.3681    0.6845    0.9048    0.9980
         0    0.4818    0.8443    0.9980    0.9048
         0    0.5878    0.9511    0.9511    0.5878

  Columns 6 through 10

    0.9511    0.7705    0.4818    0.1253   -0.2487
    0.5878    0.1253   -0.3681   -0.7705   -0.9823
    0.0000   -0.5878   -0.9511   -0.9511   -0.5878

  Columns 11 through 15

   -0.5878   -0.8443   -0.9823   -0.9823   -0.8443
   -0.9511   -0.6845   -0.2487    0.2487    0.6845
   -0.0000    0.5878    0.9511    0.9511    0.5878

  Columns 16 through 20

   -0.5878   -0.2487    0.1253    0.4818    0.7705
    0.9511    0.9823    0.7705    0.3681   -0.1253
    0.0000   -0.5878   -0.9511   -0.9511   -0.5878

  Columns 21 through 25

    0.9511    0.9980    0.9048    0.6845    0.3681
   -0.5878   -0.9048   -0.9980   -0.8443   -0.4818
   -0.0000    0.5878    0.9511    0.9511    0.5878

  Columns 26 through 30

    0.0000   -0.3681   -0.6845   -0.9048   -0.9980
   -0.0000    0.4818    0.8443    0.9980    0.9048
    0.0000   -0.5878   -0.9511   -0.9511   -0.5878

  Columns 31 through 35

   -0.9511   -0.7705   -0.4818   -0.1253    0.2487
    0.5878    0.1253   -0.3681   -0.7705   -0.9823
   -0.0000    0.5878    0.9511    0.9511    0.5878

  Columns 36 through 40

    0.5878    0.8443    0.9823    0.9823    0.8443
   -0.9511   -0.6845   -0.2487    0.2487    0.6845
    0.0000   -0.5878   -0.9511   -0.9511   -0.5878

  Columns 41 through 45

    0.5878    0.2487   -0.1253   -0.4818   -0.7705
    0.9511    0.9823    0.7705    0.3681   -0.1253
   -0.0000    0.5878    0.9511    0.9511    0.5878

  Columns 46 through 50

   -0.9511   -0.9980   -0.9048   -0.6845   -0.3681
   -0.5878   -0.9048   -0.9980   -0.8443   -0.4818
    0.0000   -0.5878   -0.9511   -0.9511   -0.5878

% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

curves3d
x=linspace(0,2*pi,10);
y=sin(x);
bspline3(x,y)

ans = 

      form: 'B-'
     knots: [1x14 double]
     coefs: [1x10 double]
    number: 10
     order: 4
       dim: 1

B=bspline3(x,y);
B.knots

ans =

  Columns 1 through 5

         0         0         0         0    1.3963

  Columns 6 through 10

    2.0944    2.7925    3.4907    4.1888    4.8869

  Columns 11 through 14

    6.2832    6.2832    6.2832    6.2832

x

x =

  Columns 1 through 5

         0    0.6981    1.3963    2.0944    2.7925

  Columns 6 through 10

    3.4907    4.1888    4.8869    5.5851    6.2832

B.coefs

ans =

  Columns 1 through 5

         0    0.4830    1.1127    0.9390    0.3711

  Columns 6 through 10

   -0.3711   -0.9390   -1.1127   -0.4830   -0.0000

xx=linspace(0,2*pi,100);
plot(xx,bspline3(x,y,xx))
plot(xx,bspline3(x,y,xx),x,y,'*')
norm(bspline3(x,y,xx)-spline(x,y,xx),inf)

ans =

   3.3307e-16

help csapi
 CSAPI Cubic spline interpolant with not-a-knot end condition.
 
    PP  = CSAPI(X,Y)  returns the cubic spline interpolant (in ppform) to the 
    given data (X,Y) using the not-a-knot end conditions.
    The interpolant matches, at the data site X(j), the given data value
    Y(:,j), j=1:length(X). The data values may be scalars, vectors, matrices,
    or even ND-arrays. Data points with the same site are averaged.
    For interpolation to gridded data, see below.
 
    CSAPI(X,Y,XX)  is the same as FNVAL(CSAPI(X,Y),XX).
 
    For example, 
 
       values = csapi([-1:5]*(pi/2),[-1 0 1 0 -1 0 1], linspace(0,2*pi));
 
    gives a surprisingly good fine sequence of values for the sine over its
    period.
 
    It is also possible to interpolate to data values on a rectangular grid,
    as follows:
 
    PP = CSAPI({X1, ...,Xm},Y)  returns the m-cubic spline interpolant (in
    ppform) that matches the data value Y(:,j1,...,jm) at the data site
    (X1(j1), ..., Xm(jm)), for ji=1:length(Xi) and i=1:m.
    The entries of each Xi must be distinct.  Y must have size
    [d,length(X1),...,.length(Xm)], with  d  a vector of natural numbers, and
    with an empty  d  acceptable when the function is to be scalar-valued.
 
    CSAPI({X1, ...,Xm},Y,XX)  is the same as FNVAL(CSAPI({X1,...,Xm},Y),XX).
 
    For example, the statements
 
       x = -1:.2:1; y=-1:.25:1; [xx, yy] = ndgrid(x,y); 
       z = sin(10*(xx.^2+yy.^2)); pp = csapi({x,y},z);
       fnplt(pp)
 
    produce the picture of an interpolant to a bivariate function. 
    Use of MESHGRID instead of NDGRID here would produce an error.
 
    See also <a href="matlab:help csape">csape</a>, <a href="matlab:help spapi">spapi</a>, <a href="matlab:help spline">spline</a>, <a href="matlab:help ndgrid">ndgrid</a>.

    Reference page in Help browser
       <a href="matlab:doc csapi">doc csapi</a>

help spapi
 SPAPI Spline interpolation.
 
    SPAPI(KNOTS,X,Y)      (with both KNOTS and X vectors)
    returns the spline  f  (if any) of order
          k := length(KNOTS) - length(X)
    with knot sequence KNOTS for which Y(:,j) equals f(X(j)) , all j.
    This is taken in the osculatory sense in case some data sites are
    repeated, i.e., in the sense that  D^m(j) f(X(j)) = Y(:,j)  with
    m(j) := #{ i<j : X(i) = X(j) }, and  D^m(j) f  the m(j)-th derivative of f.
    If you don't want this interpretation, call SPAPI with an additional
    (fourth) input argument in which case the average of all data values
    at the same data site is matched at that site.
    The data values, Y(:,j), may be scalars, vectors, matrices, or even
    ND-arrays.
 
    SPAPI(K,X,Y)          (with K a positive integer)
    is the same as  SPAPI(APTKNT(sort(X),K),X,Y), i.e., it uses APTKNT to
    obtain from K and X a suitable knot sequence.
 
    SPAPI({KNOTS1,...,KNOTSm},{X1,...,Xm},Y)
    returns in SP the m-variate tensor-product spline of coordinate order
       ki := length(KNOTSi) - length(Xi)
    with knot sequence KNOTSi in the i-th variable, i=1,...,m, for which
       Y(:,j1,...,jm) = f(X1(j1),...,Xm(jm)),  all j := (j1,...,jm) .
    As in the univariate case, KNOTSi may also be a positive integer, in
    which case the i-th knot sequence is obtained from Xi via APTKNT.
    Note the possibility of interpolating to d-valued data. However,
    in contrast to the univariate case, if the data to be interpolated are
    scalar-valued, then the input array Y is permitted to be m-dimensional,
    in which case
    Y(j1,...,jm) = f(X1(j1),...,Xm(jm)),  all j := (j1,...,jm) .
    Multiplicities in the sequences Xi lead to osculatory interpolation just
    as in the univariate case.
 
    For example, if the points in the vector  t  are all distinct, then
 
       sp = spapi(augknt(t,4,2),[t t],[cos(t) -sin(t)]);
 
    provides the C^1 piecewise cubic Hermite interpolant to the function
    f(x) = cos(x)  at the points in  t . If matching of slopes is only
    required at some subsequence  s  of  t  but that includes the leftmost
    and the rightmost point in  t , one would use instead
 
       sp = spapi( augknt([t s],4), [t s], [cos(t) -sin(s)] );
 
    or even just
 
       sp = spapi( 4, [t s], [cos(t) -sin(s)] );
 
    and the last works even if s fails to include the extreme points of t,
    and produces even a C^2 piecewise cubic interpolant to these Hermite data.
 
    As another example,
 
       sp = spapi( {[0 0 1 1],[0 0 1 1]}, {[0 1],[0 1]}, [0 0;0 1] );
 
    constructs the bilinear interpolant to values at the corners of a square,
    as would the statement
 
       sp = spapi({2,2},{[0 1],[0 1]},[0 0;0 1]);
 
    As another example, the statements
 
       x = -2:.5:2; y=-1:.25:1; [xx, yy] = ndgrid(x,y);
       z = exp(-(xx.^2+yy.^2));
       sp = spapi({3,4},{x,y},z);
       fnplt(sp)
 
    produce the picture of an interpolant (piecewise quadratic in x,
    piecewise cubic in y) to a bivariate function.
    Use of MESHGRID instead of NDGRID here would produce an error.
 
    As an illustration of osculatory interpolation to gridded data, here
    is complete bicubic interpolation, with the data explicitly derived
    from the bicubic polynomial  g(x,y) = x^3y^3, to make it easy for
    you to see exactly where the slopes and slopes of slopes (i.e., cross
    derivatives) must be placed in the data values supplied. Since our  g  is
    a bicubic polynomial, its interpolant, f , must be  g  itself.
    We test this.
 
       sites = {[0,1],[0,2]}; coefs = zeros(4,4); coefs(1,1) = 1;
       g = ppmak(sites,coefs);
       Dxg = fnval(fnder(g,[1,0]),sites);
       Dyg = fnval(fnder(g,[0,1]),sites);
       Dxyg = fnval(fnder(g,[1,1]),sites);
       f = spapi({4,4}, {sites{1}([1,2,1,2]),sites{2}([1,2,1,2])}, ...
                [fnval(g,sites), Dyg ; ...
                 Dxg.'         , Dxyg]);
       if any( squeeze( fnbrk(fn2fm(f,'pp'), 'c') ) - coefs )
         'something went wrong', end
 
    SPAPI(...,'noderiv') does not do osculatory interpolation but, instead, 
    merely averages all data values at the same site.
 
    See also <a href="matlab:help spaps">spaps</a>, <a href="matlab:help spap2">spap2</a>.

    Reference page in Help browser
       <a href="matlab:doc spapi">doc spapi</a>

diary off