% Several boundary-valued problems.

% Francisco J. Beron-Vera, 2006/15/05

%% u''(x) + u(x) = 0, u(0) = 1, u(1) = 2 (LINEAR)

% FINITE DIFFERENCES
% U(I) = u((I-1)*Dx), i = 1,...,N+1
% (U(I+1) - 2*U(I) + U(I-1))/Dx^2 + U(I) = 0, U(1) = 1, U(N+1) = 2
% E.g., N = 6:
%            |-2  1  0  0  0||U(2)|   |U(2)|                |U(0)=1|
%            | 1 -2  1  0  0||U(3)|   |U(3)|                |   0  |
% (1/Dx^2) * | 0  1 -2  1  0||U(4)| + |U(4)| = - (1/Dx^2) * |   0  | 
%            | 0  0  1 -2  1||U(5)|   |U(5)|                |   0  |
%            | 0  0  0  1 -2||U(6)|   |U(6)|                |U(7)=2|

clear

N = 50; Dx = 1/N; x = (0:N)*Dx;
d2dx2 = (diag(ones(N-2,1),1) - 2*eye(N-1) + diag(ones(N-2,1),-1))/Dx^2;
u = [1; (d2dx2 + eye(N-1))\[-1/Dx^2; zeros(N-3,1); -2/Dx^2]; 2];

figure(1), clf
plot(x,u,'o'), hold on
plot(x,cos(x)+(2-cos(1))/sin(1)*sin(x),'r.-'), hold off
axis tight
xlabel('$x$',...
   'FontSize',14,...
   'Interpreter', 'latex')
ylabel('$u(x)$',...
   'FontSize',14,...
   'Interpreter', 'latex')

% SHOOTING
% v''(x) + v(x) = 0, v(0) = 1, v'(0) = y, E(y) = v(1) - 2 (cf. SHOOT.M)

clear

[x u y E] = shoot(@(x,z) [z(2); -z(1)],linspace(0,1,25),1,2);

figure(2), clf
subplot(221)
plot(x,u,'ob'), hold on
plot(x,cos(x)+(2-cos(1))*sin(x)/sin(1),'r.-'), hold off
axis([0 1 0 3])
xlabel('$x$',...
   'FontSize',14,...
   'Interpreter','latex')
ylabel('$u(x)$',...
   'FontSize',14,...
   'Interpreter','latex')
subplot(222)
v1=y(end);
[y I]=sort(y); E=E(I);
plot(y,E,'b',[y(1) y(end)],[0 0],'k',[v1 v1],[-1 1],':r')
axis([-inf inf -1 1])
set(gca,'XTick',v1)
xlabel('$y$',...
   'FontSize',14,...
   'Interpreter','latex')
ylabel('$E(y)$',...
   'FontSize',14,...
   'Interpreter','latex')

%% u''(x)  - 0.01 * u(x)^2  + 9.8 = 0, u(0) = 0, u(5) = 40 (NONLINEAR)

clear

% FINITE DIFFERENCES
% U(I) = u((I-1)*Dx), I = 1,...,N+1
% (U(I+1) - 2*U(I) + U(I-1))/Dx^2 + U(I) = 0, U(1) = 0, U(N) = 40
% E.g., N = 5:
%                                |-2  1  0  0  0||U(2)|         |U(2)^2|              |U(1)=0|                 
%                                | 1 -2  1  0  0||U(3)|         |U(3)^2|              |   0  |                 
%  F(U(2),...,U(6)) = (1/Dx^2) * | 0  1 -2  1  0||U(4)| - .01 * |U(4)^2| + (1/Dx^2) * |   0  | + 9.8 = 0
%                                | 0  0  1 -2  1||U(5)|         |U(5)^2|              |   0  |                
%                                | 0  0  0  1 -2||U(6)|         |U(6)^2|              |U(7)=0|               

N = 25; Dx = 5/N; x = (0:N)*Dx;
d2dx2 = (diag(ones(N-2,1),1) - 2*eye(N-1) + diag(ones(N-2,1),-1))/Dx^2;
F = @(u) d2dx2*u - .01*u.^2 + [0; zeros(length(u)-2,1); 40/Dx^2] + 9.8;
dFdu = @(u) d2dx2 - .02*diag(u);
u = [0; newtonsys(F,dFdu,linspace(0,40,N-1)'); 40];

figure(3), clf
subplot(221)
plot(x,u,'r.-'), hold on

% SHOOTING 
% v''(x) - 0.01 * v(x)^2 + 9.8 = 0, v(0) = 0, v'(0) = y, E(y) = v(5) - 40 (cf. SHOOT.M)

[x u y E] = shoot(@(x,z) [z(2); .01*z(1)^2-9.8],linspace(0,5,25),0,40);

subplot(221)
plot(x,u,'ob'), hold off
axis([0 5 -80 80])
xlabel('$x$',...
   'FontSize',14,...
   'Interpreter','latex')
ylabel('$u(x)$',...
   'FontSize',14,...
   'Interpreter','latex')

subplot(222)
v1=y(end);
[y I]=sort(y); E=E(I);
plot(y,E,'b',[y(1) y(end)],[0 0],'k',[v1 v1],[-500 500],':r')
axis([-inf inf -500 500])
set(gca,'XTick',v1)
xlabel('$y$',...
   'FontSize',14,...
   'Interpreter','latex')
ylabel('$E(y)$',...
   'FontSize',14,...
   'Interpreter','latex')

%% u''(x) + k^2 u(x) = 0, u(0) = u(1) = 0 (EIGENVALUE PROBLEM)

% FINITE DIFFERENCES
% U(I) = u((I-1)*Dx), i = 1,...,N+1
% (U(I+1) - 2*U(I) + U(I-1))/Dx^2 + U(I) = 0, U(1) = 0, U(N+1) = 0
% E.g., N = 6:
%            |-2  1  0  0  0||U(2)|           |U(2)|
%            | 1 -2  1  0  0||U(3)|           |U(3)|
% (1/Dx^2) * | 0  1 -2  1  0||U(4)| = - k^2 * |U(4)| 
%            | 0  0  1 -2  1||U(5)|           |U(5)| 
%            | 0  0  0  1 -2||U(6)|           |U(6)|
     

clear

N = 50; Dx = 1/N; x = (0:N)*Dx;
d2dx2 = (diag(ones(N-2,1),1) - 2*eye(N-1) + diag(ones(N-2,1),-1))/Dx^2;
[u D] = eig(-d2dx2); 
[k I] = sort(diag(sqrt(D))); u = u(:,I); 
u = u./repmat(max(u),[N-1 1]); 
I = find(u(1,:) < 0); u(:,I) = -u(:,I);
u = [zeros(1,N-1); u; zeros(1,N-1)];

disp(' n       k_n         n*pi')
disp(sprintf('%2i %12.8f  %12.8f\n',[(1:4); k(1:4)'; (1:4)*pi]))

figure(4), clf
plot(x,u(:,1:4),'o-'), hold on
plot(x,sin((1:4)'*pi*x),'.-'), hold off
axis tight
xlabel('$x$',...
   'FontSize',14,...
   'Interpreter', 'latex')
ylabel('$u_n(x)$',...
   'FontSize',14,...
   'Interpreter', 'latex')

%% u''(x) + k^2 u(x) = 0, u'(0) = u'(1) = 0 (EIGENVALUE PROBLEM)

% FINITE DIFFERENCES
% U(I) = u((I-1)*Dx), i = 1,...,N+1
% (U(I+1) - 2*U(I) + U(I-1))/Dx^2 + U(I) = 0, U(1) = U(2), U(N) = U(N+1)
% E.g., N = 6:
%            |-2  1  0  0  0||U(2)|           |U(2)|
%            | 1 -2  1  0  0||U(3)|           |U(3)|
% (1/Dx^2) * | 0  1 -2  1  0||U(4)| = - k^2 * |U(4)| 
%            | 0  0  1 -2  1||U(5)|           |U(5)| 
%            | 0  0  0  1 -2||U(6)|           |U(6)|


clear

N = 50; Dx = 1/N; x = (0:N)*Dx;
d2dx2 = (diag(ones(N-2,1),1) - 2*eye(N-1) + diag(ones(N-2,1),-1))/Dx^2; 
d2dx2(1,1) = -1/Dx^2; 
d2dx2(end,end) = -1/Dx^2;
[u D] = eig(-d2dx2); 
[k I] = sort(diag(sqrt(D))); u = u(:,I); 
u = u./repmat(max(u),[N-1 1]); 
I = find(u(1,:) < 0); u(:,I) = -u(:,I);

disp(' n       k_n         n*pi')
disp(sprintf('%2i %12.8f  %12.8f\n',[(0:4); k(1:5)'; (0:4)*pi]))

figure(5), clf
plot(x(2:end-1),u(:,1:5),'o-'), hold on
plot(x,cos((0:4)'*pi*x),'.-'), hold off
axis tight
xlabel('$x$',...
   'FontSize',14,...
   'Interpreter', 'latex')
ylabel('$u_n(x)$',...
   'FontSize',14,...
   'Interpreter', 'latex'