%  Start off with function m-file  p.m generated last time, in order to 
%  illustrate minimum requirements for a function m-file:

>> type p

function [v,dv] = p(a,x)

dv = 0;
v = repmat(a(1),size(x));
for k=2:length(a)
    dv = dv.*x + v;
    v = v.*x + repmat(a(k),size(x));
end

% The first thing we notice is that this function m-file provides no help,
% i.e., it fails to start off with a comment section telling us about 
% this command, as is also clear from the response when we ask for help:

>> help p

No help comments found in p.m.

% I'll use the matlab editor to correct this:

>> edit p

%  ... then we get a better help:

>> type p

function [v,dv] = p(a,x)
%
%    [v,dv] = p(a,x)
%
%   value and derivative at x of the polynomial whose coefficients
%   are in  a .
:
dv = 0;
v = repmat(a(1),size(x));
for k=2:length(a)
    dv = dv.*x + v;
    v = v.*x + repmat(a(k),size(x));
end

>> help p

 
     [v,dv] = p(a,x)
 
    value and derivative at x of the polynomial whose coefficients
    are in  a .

%  Note the use of pointwise multiplication to provide simultaneously the 
% value of the polynomial at all the entries of the input array  x . This
% requires us to initialize  v  as an array fo the same size as  x , with
% all its entries equal to  a(1). This is done with the command  repmat.

% Note that  repmat(B,m,n) gives the matrix of size m,n obtained by putting
% B at each of its entries. E.g.,

>> repmat([1 2],2,3)
ans =
   1  2  1  2  1  2
   1  2  1  2  1  2

% Note also how dv is still initialized as the scalar   0 , but it becomes
% an array of the same size as  x  after the first assignment statement in the
% for-loop.
% Note also that this for-loop is hard to avoid, since it implements a 
% recursion.
% Finally, note that the second call on repmat, in the loop, is unnecessary
% since, in matlab, adding a scalar to an array causes that scalar to be added
% to every element of that array. In other words, we can simplify to
%    v = v.*x + a(k);


%   Now I make use of this function:

>> x = linspace(-1,1);
>> [v,dv] =p([1 0 0],x);
>> plot(x,[v dv])
??? Error using ==> plot
Vectors must be the same lengths.

% What went wrong?? I generated  x  via linspace, hence it is a row-vector.
% Therefore, [v,dv] = p(a,x) returns both v and dv as corresponding row-vectors,
% i.e., the plot statement asks to plot the 200-vector [v dv] against the
% 100-vector x . There are two ways (at least), to fix this:

>> plot(x,[v; dv]')

% ... in which the matrix [v; dv]' has the two columns, v', and  dv', and each
% is plotted against the vector x .

% The other way to fix it is to plot v and dv separately against  x :

>> plot(x,v, x,dv,'r')

%  .. and this is also fine.

% Note that matlab's polyval will work properly when given an array as input:

>> a = [1 0 0];
>> polyval(a,[-1 1;0 2])
ans =
     1     1
     0     4

% This is as it should be since matlab interprets the vector [1 0 0] as 
% specifying the  polynomial  p(x) := 1*x^{3-1} + 0*x^{3-2} + 0*x^{3-3} = x^2

% Next, had a quick glance at the section on structures and cell-arrays.
% Pointed out that structures are useful for attaching information of various
% kinds to just one variable. For example, 

>> x = linspace(0,2*pi);
>> sp = spline(x,sin(x))
sp = 
      form: 'pp'
    breaks: [1x100 double]
     coefs: [99x4 double]
    pieces: 99
     order: 4
       dim: 1

% i.e., the spline command returns a certain structure with 6 fields, and these
% can be accessed by in the form  variable.field   
% E.g.,
>> sp.form
ans =
pp

% The whole structure can be handed over to other commands, e.g., 
>> v = ppval(sp,x);
% .. for evaluation of this `spline' at the entries of the array  x .
% We will learn later that  spline(x,y) returns a cubic spline function that
% matches the given data x, y.
% So, in this particular case, we would expect  ppval(sp,x) to give us back
% the numbers  sin(x)  (up to roundoff):

>> max(abs(v - sin(x)))
ans = 4.8708e-018

%  ... and we are not disappointed.

% I was asked whether  polyfit(x,y) cared whether  x and y  are of the same
% size. I tried it out:

>> x = 1:3; y = sin(x)';
>> a = polyfit(x,y,2)
??? Error using ==> polyfit
X and Y vectors must be the same size.

>> a = polyfit(x,y',2)
a = -0.4180    1.3218   -0.0624
>> a = polyfit(x',y,2)
a = -0.4180    1.3218   -0.0624


>> xx = linspace(1,3);
>> plot(xx,polyval(a,xx),x,y','or')
hold on
>> plot(xx,sin(x),'linew',2.5)
??? Error using ==> plot
Vectors must be the same lengths.

% (I made the mistake of trying to plot xx  against sin(x) )

>> plot(xx,sin(xx),'linew',2.5)

%  ... and that is ok now.

% After that, I began to discuss polynomial interpolation, stressing that
% with  x  and  y  vectors of the same size, and of length  n , and assuming
% that the entries of  x  are all different, the command
% polyfit(x,y,n-1) returns the unique polynomial of degree < n that matches
% these data.

diary off