function [y0,y1] = airy_ts(x,P)
  
% Return the Taylor series about x0=0 for two independent solutions
% of Airy's equation y" - x*y = 0:
%  
% y0(x) = sum_{p=0}^P b_p x^{3p}
% y1(x) = sum_{p=0}^P c_p x^{3p+1}  
%
% x may be a vector of values.
  
  nx = length(x);

  % Start with term p=0
  bp = 1;
  cp = 1;
  x3p = ones(1,nx);
  x3p1= x;
  y0 = bp*x3p; 
  y1 = cp*x3p1;

  % Add other terms using Airy TS recursion formulas

  x3 = x.^3;
  for p = 1:P
    bp = bp/(3*p * (3*p-1));
    cp = cp/((3*p+1) * 3*p);
    x3p = x3p.*x3;
    x3p1 = x3p1.*x3;
    y0 = y0 + bp*x3p;
    y1 = y1 + cp*x3p1;
  end
  
    
  
  
  