function prob2()
%PROB2 plots solution to homework problem #2
%    PROB2() solves the problem on domain [0,10].  The solution y(0) and 
%    derivative y'(0) at x=0 should be supplied in y0 and yp0.  
%
%    Math 527 students : Compare the series solution to the exact 
%    solution and a solution obtained using ODE45.  
%
%    See also ode45, airy, complex.
%


math527 = false;   % Set to 'true' if you are a Math 527 student

clf  % Clear figure window

% --------------------------
% User defined problem setup
% --------------------------
b = 10;    % Domain : [0,b]
y0 = 1;    % y(0)
yp0 = 0;   % y'(0)

% --------------------
% Numerical parameters
% --------------------
N = 500;
x = linspace(0,b,N)';    % Domain

% ---------------
% Series solution
% ---------------
a = generate_coeffs(y0,yp0);

check_your_coeffs(a,y0,yp0);

% 5 term series
ac = fliplr(a(1:6));  % Put coefficients in descending order
y_series = polyval(ac,x);
p(1)= plot(x,y_series,'m','linewidth',3);
lstr{1} = 'Series solution (5 terms)';
hold on;

% 10 term series
ac = fliplr(a(1:11));  % Put coefficients in descending order
y_series = polyval(ac,x);
p(2)= plot(x,y_series,'b','linewidth',3);
lstr{2} = 'Series solution (10 terms)';

% 100 term series
ac = fliplr(a(1:101));  % Put coefficients in descending order
y_series = polyval(ac,x);
p(3)= plot(x,y_series,'g','linewidth',3);
lstr{3} = 'Series solution (100 terms)';

if math527
    % ------------
    % ODE solution
    % ------------
    hold on;
    Y0 = [y0;yp0];
    opt.AbsTol = 1e-12;
    opt.RelTol = 1e-12;
    [~,yout] = ode45(@frhs,x,Y0,opt);
    y_solve = yout(:,1);
    p(4) = plot(x,y_solve,'r--','linewidth',2);
    lstr{4} = 'ODE solver solution';

    % --------------
    % Exact solution
    % --------------
    y_exact = ytrue(x,y0,yp0);
    p(5) = plot(x,y_exact,'k.','markersize',6);
    lstr{5} = 'Exact solution';
end
    
% Title, labels, legend
lh = legend(p,lstr,'location','northwest');
set(lh,'fontsize',16);
xlabel('x','fontsize',16);
ylabel('y(x)','fontsize',16);
tstr = sprintf('y'''' - y'' + xy = 0');
title(tstr,'fontsize',18);

axis([0,b,-40,70]);

shg   % Bring current figure window to the front
    
% Print PNG file to turn in.
% print -dpng prob2.png


end

% -----------------------
% Check your coefficients
% -----------------------
function check_your_coeffs(a,y0,yp0)

a0 = y0;
a1 = yp0;

aa(1) = a0;                        % a0
aa(2) = a1;                        % a1
aa(3) = 0.5*a1;                    % a2
aa(4) = 1/factorial(3)*(a1-a0);    % a3

error('TODO : Supply analytic expressions for coeffs a4 - a10');
% aa(5) = ...
% aa(6) = ...
% aa(7) = ...
% aa(8) = ...
% aa(9) = ...
% aa(10) = ...
% aa(11) = ...


check_coefficients(a,aa);

end


function check_coefficients(a,aa)
% Nothing to do here

fprintf('%2s %24s %24s %12s\n','m','Computed','Your coeffs.','error');
l = double('-')*ones(1,65);
fprintf('%s\n',l);
emax = 0;
for i = 1:length(aa)
    m = i - 1;
    e = abs(a(i) - aa(i));
    emax = max([e,emax]);
    fprintf('%2d %24.16e %24.16e %12.2e\n',m, a(i),aa(i),e);
end
if (emax < eps(1))
    fprintf('Your coefficients look good!\n');
else
    fprintf('Double check your coefficients\n');
end



end


% ----------------------------
% Generate series coefficients
% ----------------------------
function a = generate_coeffs(y0,yp1)

% Use series solution to figure out how to use y(0) and 
% y'(0). 
error('TODO : Use y(0) and y''(0) to define a0, a1');

a(1) = a0;      % a0
a(2) = a1;      % a1
a(3) = 0.5*a1;  % a2

M = 100;    % Maximum number of terms 
for i = 2:M+1
    m = i-1;
    error('TODO : Define recursion relation');
end

end

% ------------------------------
% Right hand side for ODE solver
% ------------------------------
function f = frhs(x,u)

error('TODO (Math 527) : Define rhs for ODE solver');

f = f(:);  % return column vector
end

% --------------
% Exact solution
% --------------
function y = ytrue(x,y0,yp0)

error('TODO (Math 527) : Define exact solution (use WolframAlpha)');

end