% Script to plot Frobenius series y0(x,P) and y1(x,P) truncated
% after term P.

% User-chosen parameters
  xl = 0;  % left limit of plot
  xr = 10;   % right limit of plot
  Pchoices = [1 5 9 13]; % P's to plot

  x = xl + (0:0.001:1)*(xr-xl); % Plotting points
  nx = length(x);
  nP = length(Pchoices);

  color = 'rmgbk';  
  for j = 1:nP
    P = Pchoices(j);
    [J0,Y0] = bessel_frobenius(x,P);
    plot(x,J0,[color(j) '-']) % cycle through plot colors
    Y0P(j,1:nx) = Y0;
    hold on
    legend_entry(j)= {['P=' int2str(P)]}; % j'th cell of a cell array
  end

  ylim([-2 2])
  xlabel('x')
  ylabel('y_1^P(x)')
  legend(legend_entry{1:nP},'Location','SouthEast')

  for j = 1:nP
    plot(x,Y0P(j,:),[color(j) '--']) % cycle through plot colors
  end

  title(['Truncated Frobenius Series approximants to J_0(x) (solid)' ...
	 ' and Y_0(x) (dashed)'])
  hold off