% plot_Greens_fcn_1D.m

%

% This MATLAB script file makes a plot of the

% Green's function for the 1-D diffusion equation

% with the boundary conditions u(0)=u(L)=0.

%

% Both the form of the Green's function obtained

% by eigenfunction expansion and that obtained

% by the classical method are plotted.

%

% K. Beers

% MIT ChE 7/28/2002

% First, set the problem parameters.

L = 1;  % length of domain

N = 100;  % truncation order of eigenfunction expansion

zeta = 0.25;  % value of zeta for which to calculate Green's function

% Create a vector of x-values

x = linspace(0,L,100);

% Calculate the Green's function at each point using the

% classical formula.

g_class = 0*x;

for i=1:length(x)

    if(x(i) <= zeta)

        g_class(i) = -(L-zeta)*x(i)/L;

    else

        g_class(i) = -zeta*(L-x(i))/L;

    end

end

% Next, calculate the Green's formula from the eigenfunction

% expansion.

g_eigen = 0*x;

for m=1:N

    var1 = 2*L*sin(m*pi*zeta/L)/(m^2 * pi^2);

    for i=1:length(x)

        g_eigen(i) = g_eigen(i) - ...

            var1*sin(m*pi*x(i)/L);

    end

end

% Now, plot the two Green's functions.

figure;

subplot(2,1,1);

plot(x,g_class);

xlabel('x');

phrase_y = ['g(x,\zeta=', num2str(zeta), ')'];

ylabel(phrase_y);

title('Green''s function, classical (upper), eigenfunction (lower)');

subplot(2,1,2);

plot(x,g_eigen);

xlabel('x');

ylabel(phrase_y);

%

phrase_DE = 'd^2g/dx^2 = \delta(x-\zeta)';

gtext(phrase_DE);

phrase_BC = 'g(0) = 0, g(L) = 0';

gtext(phrase_BC);

%

phrase_eigen = ['expansion terminated at m = ', int2str(N)];

gtext(phrase_eigen);