function [x,t,u] = MOL_Euler(d0,d1,f,G,T,h,m)

% Use the method of lines to solve
%    u_t = u_xx + G(x,t), 0 < x < 1, 0 < t < T
% with boundary conditions
%    u(0,t) = d0(t), u(1,t) = d1(t)
% and initial condition
%    u(x,0) = f(x).
% Use Euler's method to solve the system of ODEs.
% For the discretization, use a spatial stepsize of
% delta=1/m and a time step of h. 
%
% For numerical stability, use a time step of 
% h = 1/(2*m^2) or smaller.

x = linspace(0,1,m+1)'; delta = 1/m; delta_sqr = delta^2;
t = (0:h:T)'; N = length(t);

% Initialize u.
u = zeros(m+1,N);
u(:,1) = f(x);
u(1,:) = d0(t); u(m+1,:) = d1(t);

% Solve for u using Euler's method.
for n=1:N-1
    g = G(x(2:m),t(n));
    u(2:m,n+1) = u(2:m,n) + (h/delta_sqr)*(u(1:(m-1),n) ...
                             - 2*u(2:m,n) + u(3:(m+1),n)) + h*g;
end
u = u';
end % MOL_Euler