% This script compares FD and PS solution of the 
% linear 2-point BVPs of the form
% y''(t) = exp(4*t)
% with BCs y(-1) = 0, y(1) = 0
% PS code based on Trefethen: Spectral Methods in Matlab
clear all
clc
close all
warning off
scrnsz = get(0, 'Screensize');

u = inline('exp(4*t)');
v = inline('zeros(size(t))');
w = inline('zeros(size(t))');

a = -1;
b = 1;
alpha = 0;
beta = 0;

%%%% exact solution
tt = linspace(a,b,101);
exact = ( exp(4*tt) - sinh(4)*tt - cosh(4) )/16;

%%%% FD solution
m = 99;
[t,y] = bvpsolve(u,v,w,a,b,alpha,beta,m);
figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
hold on
plot(t,y,'r')
plot(tt,exact,'g')
title([num2str(m+2) ' point FD: max err = ' num2str(norm(y'-exact,inf))],'fontsize',12)
hold off
pause

%%%% PS solution
N = 16;
[t,y] = psbvp(u,N);
% for nonzero alpha, beta use
% [t,y] = psbvp2(u,N,alpha,beta);

figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
plot(t,y,'.','markersize',16)
yy = polyval(polyfit(t,y,N),tt);    % interpolate grid data
line(tt,yy)
title([num2str(N+1) ' point PS: max err = ' num2str(norm(yy-exact,inf))],'fontsize',12)
pause

%%%% an accurate FD solution
m = 99999;
tt = linspace(a,b,m+2);
exact = ( exp(4*tt) - sinh(4)*tt - cosh(4) )/16;

[t,y] = bvpsolve(u,v,w,a,b,alpha,beta,m);
figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
hold on
plot(t,y,'r')
plot(tt,exact,'g')
title([num2str(m+2) ' point FD: max err = ' num2str(norm(y'-exact,inf))],'fontsize',12)
hold off

